Resonant Processes in a Frozen Gas 



O 
O 
(N 

u 

I Jeremy Shane Frasier 

. Charlottesville, Virginia 



C/3 



B.S. University of Virginia, 1996 



^ ' University of Virginia in Candidacy for the Degree of Doctor of 

Philosophy 



A Dissertation presented to the Graduate Faculty of the 



in 

. Department of Physics 



B 

' University of Virginia 

Ch ! May, 2001 

O 
O 



X 



© Copyright by 
Jeremy Shane Frasier 
All Rights Reserved 
May, 2001 



Processes in a Frozen Gas 



Acknowledgements 

First I would like to thank my advisor, Vittorio Celli, for agreeing 
to take me on as a student and for finding an interesting and novel 
problem on which to work. 

Thanks to Michael Fowler, we were lucky enough to have Tom 
Blum work with us during the summer of 1998. We made a great 
deal of progress and submitted our first publication on this work 
that summer, in large part because of Tom's diligent efforts. He 
also kept me from becoming lazy while my advisor was away. 

Tom Gallagher first brought the results of his experiments to our 
attention, and he was always willing to discuss this work with us. 
He also provided me with the opportunity to travel to the Univer- 
sity of Paris and work with Vladimir Akulin during the spring of 
1999. 

Ziad Maassarani, Tim Newman, and Mike Robinson provided many 
useful discussions. Bryan Wright was always willing to lend his 
considerable computer expertise. Bobby Anderson and Jay Lowell 
kindly provided the raw data for the figures that appear in their 
theses. Jay also donated a postscript file he created of the en- 
ergy levels of rubidium, which saved me the enormous trouble of 
recreating that figure from scratch. 

Acknowledgement is made to the Thomas F. and Kate Miller Jef- 
fress Memorial Trust for the support of this research. I was sup- 
ported by a Dissertation Year Fellowship from the Graduate School 
of Arts and Sciences during the fall and spring semesters of the 
1999-2000 school year. During the 1996-1999 school years I was 
supported by a Presidential Fellowship and a Hugh P. Kelly Fel- 
lowship from the Department of Physics. 



Processes in a Frozen Gas 



Abstract 

In this thesis are presented analytical results and numerical simula- 
tions concerning resonant processes in a frozen gas. This work was 
begun with the intent of modeling experiments involving energy 
transfer due to dipole-dipole interactions between ultracold Ryd- 
berg atoms, and the subsequent propagation of dipolar excitons. 
However, the results and the techniques developed to obtain them 
are more generally applicable to a wide variety of problems. 

We start by studying a simplified model that considers only the en- 
ergy transfer process in an infinite medium of randomly distributed 
atoms. This model is very instructive, because it gives a simplified 
arena in which to develop the techniques that are used again and 
again throughout the thesis. Exact analytical results are obtained 
for the time evolution of energy transfer, its dependence on the de- 
tuning from resonance, and the localization distance over which the 
transfer process is effective. The result for the initial rate of energy 
transfer turns out to have general validity. This soluble model also 
provides an excellent means for testing the numerical simulations 
that are carried out for a finite sample. 

When we allow for the propagation of dipolar excitons, it is no 
longer possible to obtain exact analytical solutions, but an approx- 
imation has been found that agrees quite well with the numerical 
simulations for the amplitude of the overall process. Simulations 
are necessary to find the sample-averaged probability of the overall 
process, and also to take into account the dependence of the results 
on the sample shape. 

In the experiments we wish to model, the atomic states are de- 
generate with respect to the sign of the projection of total angular 

momentum in the direction of the applied electric field. We show 
that, in one case of interest, this degeneracy simply leads to a 
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rescaling of the dipole-dipole interaction in the final formulas for 
the simple model. 
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Chapter 1 

Introduction 



Frozen gases are a new and in many ways ideal laboratory to test our understanding of quantum 
theory in a complex system. With present technology one can manipulate and detect electronic pro- 
cesses with an extraordinary selectivity and precision. Furthermore, the translational temperature 
of the gas can be lowered to the point where it can be ignored when discussing electronic processes. 
Frozen Rydberg gases have the added advantages that their states are well understood and many 
processes occur on a microsecond time scale, easily allowing for time-resolved spectroscopy. 

Pioneering experiments on resonant processes in these gases have been carried out by Anderson et 
al. 1^, ^ and Mourachko et al. 44 . The present work was motivated by the desire to understand 
those experiments and the subsequent work of Lowell et al. I ^Of highlighting the dynamic aspects 
of these resonant, many particle systems. The measurements of Anderson et al. and Lowell et al. 
are performed on mixtures of cold (about 150 /xK) ^^Rb atoms initially prepared in the 25si/2 and 
33sx/2 states, henceforth to be called the s and s' states respectively. There are initially N atoms in 
the s state and N' atoms in the s' state. A dipole matrix element ^ connects the state s to the state 
p, and another dipole matrix element ^' connects the state s' to the state p', where p refers to the 
24pi/2 state andp' refers to the 34p3/2 state. Dipole-dipole interactions cause ss' pp' interactions 
to occur, and then sp ps and s'p' — > p's' flippings become possible.]^ The number of p' states 
is measured after a given time has elapsed, up to about 5 /is. As a function of time, the number 
of p' states rises rapidly at first and then slowly saturates. The energy of a pp' pair can be swept 

^Although the s'p' — » p's' process is present, it is smaller than the sp ps process by a factor of 16 in this case 
and can therefore be neglected to a first approximation, except when N <^ N' . 
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through resonance with an ss' pair, where eg' + £5 — ep' + Sp, via the Stark effect by the application 
of an electric field.^ At any one time the resonance lineshape can be obtained as a function of the 
detuning A — epi + ep — es> — i which is a known function of the applied electric field. 



Figure 1.1 shows the s, s', p, and p' energies as a function of the applied electric field. This 
figure is a slightly modified version of Fig. 3.1 on page 54 of Ref. Because the p' state is in the 
j — 3/2 state, there are actually two resonances. One corresponds to the two \mj \ — 1/2 states, and 
the other corresponds to the two \mj\ =3/2 states. To a first approximation the two resonances 
can be treated separately. 



Figure 1.2 is a schematic of the ss' pp' and sp ps processes. In this illustration the circles 
represent atoms in an s state and the ovals represent atoms in a p state. In addition, the black 
shapes represent primed states, while the grey shapes represent unprimed states. The first panel 
shows a single atom in the s' state surrounded by atoms in the s state. In going from the first to 
the second panel an ss' pair has made the ss' — > pp' transition via an interaction potential V . In 
the third panel an sp pair has made the sp — > ps transition via an interaction potential U . In the 
final panel the pp' pair has made the ss' — > pp' transition, and the system has returned to the initial 
state of the first panel. 

Several features of the experimental data present a challenge to theorists and are investigated in 
this thesis. The first feature we wish to understand is the rapid rise followed by slow approach to 



saturation of the signal. This behavior can be seen in Fig. 1.3 , which is a plot of the experimental data 
taken from Fig. 3.8 on page 71 of Ref. with all the data points shown. Each individual interaction 
between pairs of atoms leads to a coherent oscillatory behavior, but as we will see in Chapter ^, 
averaging the ss' — > pp' interaction over the random positions of the atoms greatly smooths out 
the signal. The effective incoherence brought about by the sp — > ps process (the "walking away" 
discussed by Mourachko et al. |45: ij]) completes the smoothing out of the on- resonance signal but 
has less effect on the off-resonance signal. 

We also want to understand the width of the lineshape in the detuning A as a function of time. 



Fig. 1.4 shows a plot of the experimental data for the lineshape as a function of detuning for an 
interaction time of 1.78 /iS. These data are taken from Fig. 3.9 on page 74 of Ref. |4^. We note the 
presence of two resonances, which as was stated before correspond to the \mj\ — 1/2 and \mj \ —3/2 
states of the p' atom. The behavior of the collisional resonance width as a function of temperature 

•^For a discussion of the Stark effect in hydrogenic atoms, sec pages 228-241 of Ref. |^ or pages 164-171 of Ref. Jis| ] . 
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Figure 1.1: The s, s' , p, and p' energy levels of Rb as a function of electric field. 
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Figure 1.2: The circles represent atoms in an s state and the ovals represent atoms in a p state. The black shapes 
arc in primed states, while the grey shapes are in unprimed states. The first panel shows a single atom in the s' 
state surrounded by atoms in the s state. In the second panel an ss' pair has made the ss' pp' transition via an 
interaction potential V. In the third panel an sp pair has made the sp — > ps transition via an interaction potential 
U. In the final panel the pp' pair has made the ss' — > pp' transition, and the system is again in its initial state. 
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Figure 1.3: A repreaentative plot of the experimental signal as a function of time. The data are taken from Fig. 3.8 
on page 71 of Ref. pcj| with all the data points shown. 
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for hot (i.e., roughly room temperature) Rydberg gases is known to be 

Ai/(xT3/4, (1.1) 



as can be found in Ref. |21 or on page 293 of Ref. |2^. This result agrees quite well with experi- 
ments on such gases,^ and is obtained by considering only two-body collisions. Simply plugging the 
temperature 150 /iK into this result gives a width that is smaller than the room temperature width 
by a factor of about 2 x 10"*. The experimental results, however, show a reduction in the width 
of only a factor of ten. Thus we can conclude that the resonance effects of the frozen gas are not 
determined by two-body interactions alone. 

It was surmised by Anderson et al. [|| that the linewidth must be of the order of the average 
interaction energy. This interaction splits the ss'-pp' degeneracy and the subsequent migration of 
the p state to other atoms broadens the energy of this "elementary excitation" into a band, as 
appropriate to an amorphous solid. Cluster calculations were performed to illustrate this band 



formation 1 40 , but averaging over atom positions was done by order-of-magnitude arguments only. 
Mourachko et al. p5| , pointed out, in connection with their experiments on a different system, 
that the "walking away" of the p excitation from its original location should be regarded as a diffusion 
process, and dynamical equations for the resonance transition in the presence of this diffusion were 
written down ||, 

The present work builds on these earlier insights, but goes considerably further in the process 
of accurately averaging over atomic positions. It turns out that this averaging (effectively, a phase 
averaging) itself produces an exp (— ^/tI^) dependence at large times, characteristic of a diffusive 
process |l^. We discuss in detail in Chapter ^ a particular case that cleanly shows the effects of 
randomness and phase averaging. In this simple case, the averaging process can be carried out ex- 
actly. For the practically important r^^ interaction potential, one can proceed in a mathematically 
elegant way and obtain closed formulas that can be evaluated by a computer package such as Math- 
ematical or in some cases by hand. It appears to be a fortunate coincidence that the potential 
accurately describes the dipole-dipole interactions that are so prevalent in nature. Some aspects of 
this approach are easily extended to more general interactions; however, the dependence has 
special properties which, when coupled with the assumption of a random distribution of atoms, lead 

^See, for example, pages 309-310 of Ref. H. 
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to relatively simple results. 

Although this work was stimulated by experiments on Rydberg atoms, it is generally applicable 
to resonances induced by dipole-dipole interactions. These may be common in highly excited gases 
and in molecular systems. Frozen, resonant gases may even exist in interstellar clouds. 

First, in Chapter |^, we discuss a single s' atom interacting resonantly with a surrounding gas 
of s atoms, without any possibility for the s atom to interact with p atoms through an sp ps 
process. When the result is averaged over the positions of the s atoms, it corresponds to a "sparse" 
system of s' atoms, i.e. a system where N' <ti N. We do not know of an existing experiment to 
which this treatment applies, so this section can be viewed as a theoretical prediction of the outcome 
of a possible experiment. Our purpose here is also one of exposition, since this example allows us 
to introduce in the simplest context some of the mathematical techniques used throughout the rest 
of this thesis. Some plots of the results are presented, and in the final sections of this chapter we 
consider the effect of interaction potentials more complicated than just a simple l/r"^ potential. 

In Chapter |^ we continue with our discussion of the simple model of Chapter ^ by examining 
the issue of localization. We study localization in this system because we need to have some idea 
of the relative importance of close pairs of atoms, as well as some idea of the effective range of 
the interaction. Naively, it seems that the close pairs are dominant and that the range of the 
l/r^ potential is infinite because l/r^ diverges at r = and yet the integral J (fir/r^ over all space 
diverges at both limits. 

In Chapter ^ we come one step closer to the actual system used in the experiments by allowing 
the sp —^ ps process to take place. This renders the problem analytically insoluble but we consider 
several approximations, each of which allows us to make some progress. When compared with the 
results of numerical simulations, we see that the coherent potential approximation compares most 
favorably. 

In Chapter ^, we introduce an effective spin Hamiltonian which fully models the Rb system of 
Anderson et al. |^ and Lowell et al. for all values of N and N' , as well as for any strength 
of the resonant process ss' —^ pp' and of the mixing processes sp ps and s'p' — > p's'. This 
Hamiltonian corresponds to a system of two interpenetrating spin glasses. We show its equivalence 
to the results of Chapter |^ in the sparse limit (N' <^ N). We then derive an entirely equivalent but 
seemingly more cumbersome Hamiltonian to describe this system. This second derivation embodies 
a more general technique, which we demonstrate by constructing a spin Hamiltonian that describes 
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the cesium system studied by Mourachko et al. |45 



In Chapter ^ we take into account fully the complication introduced by the consideration of spin, 
and we determine exactly what the spin-dependent interatomic potentials are for the sp ps and 



ss' — > pp' processes in the experiment we are modeling ||^, 40 . 

In Chapter |^ we discuss numerical simulations of the system we are considering. We start by 
describing in detail how the simulations are carried out, then we present the results, including 
comparisons of the simulations with the experimental data. 

The limitations of the present theory and some of the many possible extensions are discussed in 
Chapter |. 



Chapter 2 



Sparse ss^ pp^ Processes — The 
Simple Model 



2.1 Introduction 



In this chapter we discuss a single s' atom interacting resonantly with a surrounding gas of s atoms, 
without any possibility for the s atom to interact with p atoms through an sp ps process. In 



Section 2.2 we write down and solve the Schrodinger equation for a given spatial configuration of 



the s atoms. Then in Section 2.3 this result is averaged over the positions of the s atoms, so that it 



corresponds to a sparse system of s' atoms where N' <ti N. In Section 2.4 we compute the averaged 



signal and the lineshape. Specifically, in Sections 2.4.1 and 2.4.2 we start from the Laplace transform 



of the signal and develop series expansions for the signal as a function of time. We also discuss the use 



of the fast Fourier transform to obtain the signal as a function of time. In Section 2.4.3 we develop 



a numerical method to compute the averaged signal and lineshape that works exceedingly well, and 

we 



in Section 2.4.4 we discuss the good and bad points of the three methods. Next, in Section 2.4.5 



present plots of the averaged signal and lineshape. In Section 2.5 we consider the effect of interaction 
potentials more complicated than just a simple potential. Finally, we summarize the results of 



this chapter in Section 2.6 
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2.2 Basic equations 

We consider one atom at the origin, initially in the state s', in interaction with a gas of N randomly 

distributed s atoms through an ss' — > pp' process. Like Mourachko et al. p5[ |, we describe the system 
by the equations 

N 

iao = Aao + ^VkCk, (2.1a) 
fe=i 

ick = 14 oo. (2 -lb) 

Here ag (A, t) is the amplitude of the state in which the atom at the origin is in state s' and all 
other atoms are in state s, while Ck (A, t) (with k running from 1 to N) is the amplitude of the state 
in which the atom at the origin is in state p' and the atom at r^. is in state p, while all the others 
remain in state s. The quantity is the interaction potential and A = ej, + ep — e(, — Cs is the 
detuning from resonance. For simplicity, we will assume for now that Vk is of the dipole-dipole form 

Vk = ^, (2.2) 



and we will discuss later in Section 2.5 what happens when the interaction potential has a different 
dependence on r or contains an angular dependence. Furthermore, the atoms are sufficiently cold 
that during the time scale of interest they move only a very small fraction of their separation, and 
therefore Vk can be taken to be independent of time. As the temperature increases one approaches 
the opposite limit, where binary collisions control the resonant process ]57|. 



The set of differential equations (2.1a) and (2.1b) can easily be solved with the initial condition 



ao (A, i = 0) = 1 and Ck (A, t = 0) = for aU k, yielding 



,(^i 2.^.e-AV^sin(VA^t/2) 



with 

N 



V'-Vlf. (2.4) 



/ . "k 

k=l 
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Eq. ( 2.1b| ) then gives 



ao (A,t) = e 



COS VA2 +4V2V2 - 



VA2 + 4V2 



sin \/A2 + 4V2V2 



(2.5) 



The experiment detects transitions to any one of the p' states (or, equivalently, to any one of the p 
states), and is thus a measure of 



^(A,t) = l-|ao(A,t)|' = ^|cfe(A,i)r 

A; 

4V2 sin^ (VA2 +4V2 t/2) 



A2 + 4V2 



(2.6) 



We see that the signal S {A,t) exhibits Rabi osciUations. 

We apply this result to a system of N' atoms, initially in state s', randomly dispersed among a 
much larger number, N, of s atoms. In this sparse limit, the experimental signal is proportional to 
N' times the sample average of (A, i), and the sample average is equivalent to an ensemble average 
over the atomic positions that are hidden in (see Eq. ( |2.§| ) below). In general, averaging is 
more easily done on the Laplace transform of the signal, which in this case is 



^(A,a) 



.^(A,*) dt 



2V2 



a 



2a 



-4V2 
a2 + A2 



a2 + A2 + 4V2 



(2.7) 



This Laplace transform can also be found as Eq. (5.9) on page 28 of Ref. |51 



2.3 Averaging over atom positions 

To compute ensemble averages we use the following result, valid in the limit TV ^ 1, for a random 
distribution of the variables r^:^ 

N 

fc=i 

^The procedure here is, in effect, that used in the derivation of the second virial coefficient. See, for example, 
page 421 of Ref. 



= ^ I d ri . . .d rN exp 
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N 



1 - - I d^r 



1 

n 



1 - g-*^^ 



N 



(2.8) 



where 51 is the volume of the gas in the trap and 



(2.9) 



In particular, for the dipolar interaction V{r) — fifi' /r , we have 



1-e 



47^3/2 



leading to 



where 



4^ TV 

3 n 



fifi 



(2.10) 



(2.11) 



(2.12) 



For the typical densities in the experiments of Anderson et al. Q and Lowell et al. po|| , v is on 
the order of MHz. Since v has the units of an energy, it must be of order pp'N/D, on dimensional 



grounds. It is still remarkable, however, that Eq. (2.12) gives simply and exactly the quantity v that 
will enter in all the averaged quantities in this section.n 



Using Eq. (2.11), we can evaluate the average of any function F(V^). We simply need to 
determine the probability distribution of the variable V^, which we will denote as P(V^). By 
definition, we have 



N 



(2.13) 



We proceed by writing the delta function as a Fourier transform over the dummy variable q. This 
leads to 

(2.14) 



k=l ^ 



dipolar sums for regular lattices are discussed in Refs. 



i nfig urations c 
|l| and 0. 
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The average can now be performed using Eq. (2.11), and we have 



(2.15) 



The remaining integral requires a httle finesse to evaluate. First, we split the integral into two 
parts, so that 

(2.16) 



and 



P,{V')^J^^^ey,p{^^qV^-^qv) 



' — OO 



^°°|exp(-.,V^ 



-iqv 



(2.17) 



Now we note that 



P(V2) = Pi{V^)+P2{V 

r^dq 

lo 2^ 



2 Re / exp (^iqV^ — \/iqvj 



(2.18) 



Changing variables to cr = ^/^iqV^ ~ e ^^l^^fqx, we see that q — ia^ jV"^ and hence 



F(V2) = 



2 1 
2 1 



Re 
Re 



I \ da G exp ( —(7 



da exp — (T^ — 



/V2 



(2.19) 



This can be rewritten as 



(2.20) 



and the remaining integral is easily looked up in tables, such as Ref. p3]. One finds that 



y dx exp cos {hx) = 



(2.21) 
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and hence we have 

1 



When the derivative is evaluated we are left with 

1 V / 



It is interesting to note that P (V^) is not Gaussian, and the central limit theorem does not ap- 
ply, because P {V) cx does not have a root mean square value. In fact, P {V) is not even 
normalizable. 

Therefore, using Eq. ( p. 23 ) we can compute the average of any function F (V^) as 



{F{V^))= d{V^) P{V^)F{V^). (2.24) 
Jo 

We can further define = w^/V^, and then by a simple change of variable the average of F can be 
written in the form 

which agrees with the result derived by an alternative method in Ref. . 

This last set of relations implies that we can average a function F(V^) over the positions of the 
interacting dipoles (atoms in our case) by replacing with jy^ and integrating over the kernel 
exp (— y^/4) I \pn. This trick is not always useful as such, because it can lead to highly oscillatory 
integrals. However, if the function F has suitable analytic properties, then the integration path of 
Eq. (I^H) can be shifted in the complex y plane in such a way that the oscillations are damped out 



and the integral is rather easy to evaluate numerically. An example of this procedure is discussed 



in detail in Section 2.4.2 
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2.4 The averaged signal and lineshape 



The average over atomic positions can be obtained in closed form for the Laplace transform of the 



signal, ( S (A, a) ). Starting from Eq. (2.7) in the form 



1 + A'^ f°° 



■/3(a2+A2+4V^) 



(2.26) 



and using Eq. (2.11), we obtain 



S'(A,a) 



1. ^ 

2a ~ 2q 



d(3 e 



/TT V IV 

■ exp 



2 aA 



erfc 



(2.27) 



where A^ = + A^. For the sake of completeness, and because it will be used later on, we note 



that a similar result can be obtained for (fio (A, a)) by applying the same method to Eq. {2.1). One 
finds 



(ao (A, a)) 



1 



a + lA 



(a + iA) 



3/2 



exp 



4a (a + iA) 



erfc 



2^Ja (a + iA) 



(2.28) 



The problem is then to invert the Laplace transform of Eq. (2.27). As is shown in Section 2.4.1, 
it is possible in this case to convert the Laplace transform to a Fourier transform, which can then 
be inverted numerically by a fast Fourier transform algorithm. A more convenient inversion method 



was presented in Rcf. ||18| and is described in more detail in Section 2.4.2. It leads to an expression 



for {S (A, t)) in terms of generalized hypergeometric functions, which can be conveniently handled 
by Mathematica. 



A third method to obtain {S{A,t)) is to apply Eq. (2.25) directly to (5 (A, a)) as given by 



Eq. (2^). The remaining integral is difficult to perform numerically, as we already noted, but this 



difficulty can be circumvented as is shown in Section 2.4.3 . The advantages and limitations of the 
three numerical methods for evaluating {S (A,t)) are compared in Section ^.4.4 . 



Plots of {S {A,t)) are presented and discussed in Section 2.4.5 . Readers who are not interested 
in mathematical methods can proceed directly to Section 2.4.5| . 
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2.4.1 The Laplace and Fourier transforms 

We see from Eq. (7.1.6) on page 297 of Ref. § that 



e erfc (z) = e — > — 

^ n— 



2n+l 



Now we use the fact that 



and we find 



E 



■^r(n+l) r(2n + 2) 



E 



4"r(n + l).,2„+i 



r(2z) 



1 ^2. A 1 



(2.29) 



(2.30) 



Therefore Eq. ( 2.27 ) gives 



erfc (z) 



Er(^ + i) Er(„_, 



(2.31) 



5 (A, a) 



2 p ^n±2N I 



(2.32) 



This equation shows that (^S (A, a)j can be analytically continued to complex a with Re a > 0, 
although it has a branch cut ending at essential singularities located at a = iA and a = —iA. In 
particular, the Fourier transform is obtained by putting a — uo/i. We have obtained {S {A^t)) by 
inverting this Fourier transform numerically using a fast Fourier transform routine. However, the 
following method of inverting the Laplace transform is much more convenient, for reasons that are 



explained in Section 2.4.4 



2.4.2 Inverting the Laplace transform 

According to Eq. (2.66) of Rcf. |51|, the inverse Laplace transform of 



Q,2p (|q,2 _|_ g2'^ 



(2.33) 
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IS 

^2^+2.-1 / 1 ^2^2 



— -1F2 I V. u + 1^,11 + + -: , (2.34) 

r(2^ + 2j/)^ ^ V 2' 4 y' ^ ^ 

where the function pFq (ai, 02, . . . , ap; 61, 62, • • ■ , i'q! -z) is a generahzed hypergeometric function, and 
is defined on pages 750-751 of Ref. to be 

F(aa a-bb b ■ - HLi T (6.) f> Ut^^ {a, + k) 

(a, a2, . . . , a„ 6, 62, . . . , b„ z) - ^^^-^ D^.., r (6„ + ^) T (fc + 1) " ^^'^^^ 

Hence we find that the inverse Laplace transform of ( 5* (A, a) \ is 



(2.36) 



(^(A,t)) = — 2^(-l) -^^^^_^iF2(^^,^,^,-^ 



Expanding 1F2 using Eq. (I^H) yields 



(5(A,t)> - 2^^A.r(^)r(n + 2) r(^) 



E 



p ( »+2r»+2 -) p ( »+2r»+3 ^ T (m + 1) 



m=0 



2 ,iio r(n + 2)r(i^) 

r(^i±fi±i) / A2t2\"' 

p ^ »+2m+2 ^ p ^ n+2m+3 ^ T (m + 1) V ^ 

^ OO OO 111 



2 ^ ^ ^ ^ 2 r (n + 2) 

1 / A2<2 



p ( »+2m+2 ^ r (to + 1) V 4 



oo oo 



^ v-v- (-l)"+'"^"A''"t"+2"^ 

^,^^g4™r(n + l)r(TO + l)r(S±^!i±2)(n + 2TO+l)' ^' ^ 



which can also be obtained directly by expanding Eq. ( 2.27] ) in 1 j a and then taking the inverse 
Laplace transform term by term. 



For some purposes, it is convenient to rearrange Eq. ( ^.37 ) in the form of a series in powers of 



the detuning squared. To achieve this, we use the second line of this equation and find 

°o / a2j-2\™ 1 
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where 



n=0 



r(n + 2)r(2i±i)r(2±M)r( ^ 



r (n + I) r (n + m + ^) 



n=0 



r (2n + 2) r (n + i) r (n + m + 1) r (n + m + I) 



, — oo 



2n+l 



r(n + 2)r(n + m+ 1) 



n=0 



r (2n + 3) r (n + 1) r (n + m + I) r (n + m + 2) ■ 



(2.; 



Again using Eq. ( 2.30 ), we see that 

n=0 



r(n + |)r(n + m+i] 



9 4n+3/2 \" ^ I 



r (n + 1) r (n + I) r (n + i) r (n + m + 1) r (n + m + I) 
r(n + 2)r(n + TO + l) 



n=0 
°° / ,,2+2 



r (ti + I) r (n + 2) r (n + 1) r (n + m + I) r (n + m + 2) 
r (n + rn + 



~ 4^^„\^~y r(n+l)r(n+i)r(n + m+l)r(n + m+|) 



2 4 ^ V 4 



r (n + m + 1) 



r(n + |)r(n + 1) r (n + m + I) r (n + rn + 2) 



^ r(m + i) /II 3 i;2i2 

4 r(i)r(TO+i)r(m+f) ' ' V 2'2' 2' 4 



2j.2 
-V t 



r(m + 1) 



3 3 

, , , 1 ( TO + 1; — , m H — ,TO + 2, 

r(§)r(TO + |)r(m + 2)' H '2' 2' ' 4 



1 



11 , 3 w2i2 

I"-/ IN ,1-^3 toH — ;— ,to+1,toH — ; — 

4 (TO+i)r(TO+l) ^ V 2' 2' ' 2' 4 



2.2 



1 



1F3 m + 1; -, TO + -, TO + 2; ■ 



(TO + i)r(TO + |)'^n"' ■ ^'2''" ' 2' ■ 4 



(2.- 



Thus the expression for the signal becomes 



(^(A,i)) = 



m— 



r(TO + 1) 



1 



^,11 3 ^2^2 

-vt-. T\ -, T TO H — : — , TO + 1, TO H — : — 

4 (rri + i) r (m + 1) H 2' 2' ' 2' 4 



4 (TO + i)r(TO 



,3 3 „ t;2^2 
-1-F3 \ m+l;-,m + -,to + 2; — 



A2i2 
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-Vt- 



1 



4 - (m+l)[T{m + l)Y 
1 



1 1 



1F3 m + -; -,rn + l,rn + - 



3 i;^t^ 



2' 2 



2' 4 



, -1 F-j m + 1: — , m H — ,m + 2: 

r(m + |)r(TO + 2)' H '2' 2' '4 



(2.41) 



If we consider the signal at resonance, then the expression in Eq. (2.41) simphfies considerably. 
We note that only the m = term contributes when A = 0. Thus we have 



(^(A = 0,<)) = 



1^ 1 
^vi-. 

4 i[r(i)] 



2 

4 r(|)r(2) 



2 1^3 \ 2^ 2'-^' 2' 4 



1-F3 



/TT 



-vtoF2 



3 vH^ 
2' 4 



1 3 3 „ «^ 

'2'2' ' 4 



1,2^/^ 3 3 „ w2t2 
2-^^SF3 l;2,2,2; — 



(2.42) 



We further note that 



1-^3 



3 3 „ vH^ 



i..Vr(|)r(||r(2,. 



r(n + 1) 



r (n + 1) r {n + 1) r (n + 1) r {n + 2) 



IE 



1 



2,2 \ "+i 



2^or(« + i)r(« + i)r(n + 2) V 4 
ir(- + ^)r(n + i)r(n + i) 4 



TT 



1 



2r(i)r 



1 1 vH 



2j-2 



0^2 I ; 77, 77 

'2'2' 4 



^2^ V2 

1 /II v'^t'^\ 1 
T: 0-F2 ; 77, 77 



2' 2' 4 



(2.43) 



Thus the signal at resonance can be written as 



1 1 

2^2 



" "°'^2''2'2' 4 



1 1 vH' 



(2.44) 



which agrees perfectly with Eq. (19) of Ref. [nsl 
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2.4.3 A numerical method using Eq. (|2.25| ) 

We have from Eq. ( ^.25 ) that the averaged signal can be written as 



{S (A, t)) ^ dy^ exp (^-^) S (A, 



(2.45) 



Then, using Eq. (2_^), we have 



(5(A,t)) = — / dyexp(-i^^^^ 



^2 sin2 



A2+4^< 



A2 +44 



(2.46) 



We can rewrite this as 



{S{\t)) 



'T^ JO 

2 



dy exp 



4 y A2 + 4i4 



1 — COS 



'A2 + 4^i 

y J 



/n 



'^°^Pr4;A2 + 4- 



2 f"" 
— Re / dy 



1^ exp ( + zt Ja2 + 4^ 



(2.47) 



If A = 0, then we have 



(5(A = 0,i)) = --/(A = 0,i), 



(2.48) 



where 



I{A = 0,t) 



Re / dy exp 



' 2it- 
V 



(2.49) 



We can change the path of integration from the positive real axis to a path where we start from the 
origin and proceed along the negative imaginary axis to the point y — —i, then proceed to the point 
y — oo ~ i along the line Imy = — i. Then we have 



/(A = 0,t) 



1 



1 



Re 



1 , f^'^ 

— dz exp 2t— 

i Jo V 4 z 



20? JO 



dw exp 



{w 



2it- 



(2.50) 
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Clearly the integral over z is purely real, and so 



(5(A = 0,t)) = i--^Re / dwexp 
^ Jo 



{w — lY 



2ivt 



(2.51) 



Mathematica is quite adept at quickly evaluating this integral for any value of vt. 

The case where A is only slightly more complicated. We have from Eq. ( ^.47 ) that 



{S (A, t)) = ^^ exp (^^) erfc ( J) - / (A, t) , 



(2.52) 



where 



/A,t) = — Re/ dy—^^ exp I + ivtJ — + — ] 



(2.53) 



The path we used in the resonant case will still work well for small A/v, but will fail for A ^ 2i;. 
This is clear because the radicand ^ + -p- will change sign as y goes from zero to — i. By analogy 
with the resonant case and through a bit of trial and error, however, one finds that for A ^ 2u 
the path where we start from the origin and proceed along the positive imaginary axis to the point 
y — —-^i then proceed to the point y ~ oo ~ ^ along the line Imj/ — ^ works well. Using this 
path we have 



J(A,t) = -—Re 



1 I 1 /A2 4 



2 f°° 
—= Re / dw 

\ IT 



1 



exp 



A2 

v2 1 „2 



\2 nr^ : 



2AJ 



ivt. 



2A; 



-1 



(2.54) 



It is again the case that the z integral gives zero, and we have that 



{S{A,t)) = 



— A^xp^jerfc^-j 



2 f°° 
—= Re / dw ■ 

A ^0 



w ■ 



2A; 



exp 



2Ay 



A2 



2\/ 



(2.55) 
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Mathematica is again fairly adept at evaluating this expression numerically for any values of vt 
and A/w, although the computation time increases with increasing A/ v. In practice, we have used 



Eq. ( |2.5lD for A < and Eq. ( p.55| ) for A > t;. 

Note that, using this method, it is possible to obtain a plot of the resonance width versus time 
by solving numerically for a value of A such that 

{S{A,t)) = ^{S{0,t)), (2.56) 

at many values of t. 

We can also use this method to obtain asymptotic expressions for the signal {S (A, t)). Consider 



the resonant case. We see from Eqs. (2.48) and ( p.49| ) that we want to obtain an asymptotic 
expression for the integral 

/ (A = 0, = -^Re / dy exp [<j> (y)] , (2.57) 
W'^ Jo 



where 



, , "v^ 2ivt 

'{y)^-^ + • 2.58 

4 y 



We can simply apply the method of steepest descent to this integral, as discussed on pages 82-90 of 
Ref. II. We see that 

y 2ivt 
2 



^'(2/) = -^-9^, (2.59) 



and 



r{y) = -U^. (2.60) 

2 y-^ 



Thus we find that 0' (y) vanishes when y^ = —4ivt. To ensure convergence, we want to pick the 
root closest to the original integration path in the lower half plane, and so we choose the root 



yo = exp(-^)(4^;^)^/^ (2.61) 



We then have 



0(y)«0(yo) + ^</'"(2/o)(y-2;o). (2.62) 
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Proceeding carefully we find that 



(2.63) 



and 



r (yo) = -6. 



Finally there is only a Gaussian integral to do, and one easily finds that 



(5(A = 0,i))«l-i=cxp 



3 , ^ 



2/3 



COS 



3^/3 



(vt) 



2/3 



(2.64) 



(2.65) 



for large t. 

The same procedure can be applied when A 7^ 0, although this situation is somewhat more 
complicated because in order to determine the position of the saddle point one must first solve a 
cubic equation. This is best done numerically, and one finds that the result is qualitatively the same 



as Eq. (2.65) 



2.4.4 Comparison of the three methods 

Each of the three methods we have discussed has its drawbacks. 



In order for the fast Fourier transform method of Section 2.4.1 to be practical, we must first 
subtract off certain parts of the Fourier transform that are not well behaved. For instance, we must 
subtract off the part of the Fourier transform that goes as l/uj at small uj. Because we must sample 
points over a wide range of values of to in Fourier space in order to obtain an accurate picture of the 
signal in t space, however, we must also ensure that the Fourier transform falls off suitably fast as 
Lu increases. We actually subtract off the quantity 



exp(j^)erfc(^) 
2 Auj + i exp {-^) erfc (^) i 



(2.66) 



so that the coefficient of l/cj^ at large u! also vanishes. The parts that are subtracted off must have 
a known inverse transform, which at the end is added to the remaining part that is transformed by 
a fast Fourier transform algorithm. Finding a suitable subtraction is a rather cumbersome process. 
In addition, because the fast Fourier transform method is a numerical method (as opposed to an 
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analytical method), we do not obtain an explicit form for {S (A, t)). 

The series obtained using the method of Section 2.4.2| , on the other hand, do give an explicit 



form for (S'(A,t)). Unfortunately, however, they have the very serious limitation that they do not 
always converge for larger values of t unless A is small. 



The method discussed in Section 2.4.3 shares a drawback with the fast Fourier transform method, 
namely that unless the integral can be done analytically, an explicit form for {S {A,t)) cannot be 
obtained. This method is, however, superior to the fast Fourier transform method because it does 
not require the cumbersome subtractions. 

As a general rule, the series method of Section 2.4.2| is best when it converges, and the integral 



method of Section 2.4.3 is the best otherwise. 



2.4.5 Plots and discussion 



Plots of {S (A, t)) for several values of A are shown in Fig. 2.1. Notice that the initial slope of the 
signal is independent of A. This feature translates into resonance widths that vary as 1/t for small 
t - the so-called transform broadening discussed by Thomson et al. ^ . With v on the order of 
MHz, this initial rise occurs in a fraction of a /is. We will see later that when the sp — > ps process 
is included, the initial rise is independent of both the detuning A and the strength of the sp — > ps 
interaction. Another point of interest is that the averaging is much more effective at smoothing out 
the oscillations for A = than it is for A 7^ 0. 



We can construct a measure of the width from the first two terms in the A-expansion, Eq. (2.41). 

If 

(5 (A, t)) = 5o {t) + A^^i {t) + O (A^) , (2.67) 

and we define 



then 2w would be the FWHM if the lineshape were Lorcntzian. We note that the functions 6*0 (i) and 
Si (t) are easily obtained in closed form, since they are essentially just the m = and m = 1 terms 
in the sum of Eq. ( p. 41 ), respectively. The quantity w/v is plotted in Fig. 2.2 as the dashed curve. 



For small vt, we have explicitly w = yYljt. While the 1/t behavior follows from the transform 
broadening argument, or even simply from dimensional analysis, the coeflacient vT2 is a prediction 
of the detailed theory. 
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Figure 2Jj^_Averaged signal S (A,t) for (from top to bottom) A/v = 0, 1, 2, 4, and 8, calculated using the method of 
Section 2.4.3 for the sparse limit. The ss' — > pp' process is averaged, and the sp — > ps process is neglected. The Rabi 
oscillations are not completely washed out by the averaging process and are more pronounced off resonance. 
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The solid curve in Fig. 2.2 is the exact half width, as determined using Eq. ( p^ ) of Section |X|. 
We see that the width of Eq. ( p. 68 ) is a good approximation to the exact result only at small times. 



The lineshape at saturation (t — > oo) can be obtained directly from the first line of Eq. (2.47) 
For large times the cosine term averages to zero, giving 



(5(A,i->oo)> 



-^^expf^lerfc 



(2.69) 



which can also be extracted from the small a behavior of Eq. (2.27). This lineshape is plotted in 
Fig. The FWHM is approximately 4.6w. Also plotted is a Lorentzian lineshape with the same 
height and FWHM. Notice that this lineshape is sharper than the Lorentzian for small A and falls 
off more slowly for large A. 



2.5 Averaging more complicated interaction potentials 



In Section 2.3 we averaged the quantity e~^'^ , where = X^fc^i "^^z? — this 

section we consider the problem of averaging a more complicated interaction potential. Here we 
perform the averaging procedure for a potential that varies as 







(2.70) 



where C is a constant and r^, 9k, and (j)k are the spherical coordinates corresponding to r^. We have 
seen that we can evaluate the average of the signal if we can compute the average of the quantity 
exp (— /3V^). A careful look at Eq. ( ^.8[ ) shows that the steps do not make any assumptions on how 
V (r/j) depends on r^, so this equation still holds and we have 



exp 



N 



d^r 



1 - e-''^ W 



(2.71) 



Thus we are left with the task of evaluating the integral 



/ = 



dr r 















/ 


b I dO sin e 


1 — exp ^ 


JO 


Jo 





(2.72) 
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Figure 2.2: The Hash ed lin e is the w idth in. as given by Eq. (2.68), versus time. The solid line is the width as 
determined using Eq. ( [2.56| ) of Section |2.4.3[ The "dip" seen in bo th c urves is associated with the stronger and higher 
frequency oscillations scon at larger values of A, as shown in Fig. 2.1 
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Figure 2.3: The solixLcurve is the saturation (t — > oo) hneshape of Eq. (2.69). It is not a Lorentzian but a Lorentzian 
averaged as in Eq. (2.25). The dashed curve is a Lorentzian with the same height and FWHM, drawn for comparison. 
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We can define the variable z = [C^f^ 0)] and then 

/•TT 

d0 / dO sin 61 [C^f {d,(f>) P]^'^^ 



dz 

3 Jo J a 

i|^|3/p^3/2p 



1 — exp 



dz 



1 — exp 



r2p/3 



27r /"TT 

Jo 



(2.73) 



The remaining integrals over z and the angular variables can now be looked up in tables or done 
numerically. 

It is very important to note that the introduction of an angular dependence into the potential 
does not affect the averaging process except to change the definition of the scaling parameter v by 
a constant factor. This means that we can, for simplicity, use a simple potential that contains no 
angular dependence to get an answer, and then take into account any angular dependence by simply 
changing the definition of v in the formulas. 



2.5.1 Interaction potential of aligned dipoles 

The induced dipoles in the magneto-optical trap become aligned when the applied electric field is 
sufficiently large that the Stark effect has passed from the quadratic to the linear regime. The aligned 
dipole case is the one considered in Refs. and |5^, for example. It is also the potential we use 
in performing the numerical simulations discussed in Chapter ^ for reasons that are explained fully 



in Section 7.2. Specifically, we have 



(2.74) 



Thus we need to evaluate the integral 



1 f /""^ 

■^lJ-iP2\/Ti \j dz 



1-exp --TT 



2ir 



dO sin6' 3cos^6l- 1 



3 

16^2 
9V3 



3\/3 



-/UlM2 



(2.75) 



We see that the additional angular dependence changes v by a factor of 4/3\/3 « 0.77. 
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2.5.2 Effect of a different radial dependence on the initial rise 



In this section, we want to determine the effect of a different radial dependence in the interaction 
potential on the initial rise. We saw that if the potential goes as l/r"^, then the signal starts linearly 
with with a slope that is independent of the detuning A. We will see that if the interaction potential 
goes as 1/rP, then the initial rise will go as t^^P, with a coefficient that is again independent of the 
detuning. 

Consider again Eq. ( 2.26D in the form 



1 a' 



2a 



^0 



2a 



(2.76) 



We noted previously that any angular dependence in the interaction potential will only change the 
scale factor v, and will not change the functional form of the signal. Therefore we can simply consider 
a potential of the form 

r 

(2.77) 



Of course, we can incorporate the effect of an angular contribution by simply rescaling C. 
It follows, then, from Eqs. ( |2.7lD and ( |2.73| ) that 



(2.78) 



where 



43/2p 



dz 



1 — exp 



If we define ^ = /3 (a^ + A^) , then 



5(A,a)) = --,^ 



(a2 + A2) 



2\3/2p 



(2.79) 



(2.80) 



We want the behavior at small times, which corresponds to large a. Therefore we can replace + A^ 
with a^ and expand the second exponential, finding 
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It immediately follows that 



{S{A,t))= \ „t'^P-.... 



2r 1 



(2.82) 



Thus we see that if the interaction potential goes as l/r^, then the initial rise is independent of 
the detuning, A, and is proportional to t'^^P. As a check, we note that if p = 3, then we have 



{S{A,t)) = 



2T (2) 
—j-i^t - 



(2.83) 



where 



4^1 |C| 
•J Jo 



8^3/2 



dz 



1 — exp 



(2.84) 



and so 



{SiA,t))^y-vt 



(2.85) 



which agrees perfectly with the m = 0, n = term of Eq. (2.37) 



2.6 Summary 

In this chapter we discussed the case of a single s' atom interacting resonantly with a surrounding 
gas of s atoms, without the possibility of any sp-ps interaction. The results, when averaged over the 
positions of the s atoms, are applicable to a sparse distribution of s' atoms among the s atoms (or 
vice versa, with the appropriate changes). This case is important to study for several reasons. First, 
exact analytical results can be obtained. This not only provides an excellent means for testing the 
validity of the numerical simulations discussed in Chapter but it also provides a simpler arena for 
the development of the techniques that are used again and again in the more complicated models of 



the later chapters. Two key examples of the latter are the averaging technique of Sections 2^ and 



2.5 



and the writing of sums as generalized hypergeometrics, as was done in Section 



2.4.2 
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The averaging is non-trivial because the interaction potential itself, being proportional to 
does not have an average value. Nevertheless, an effective average potential v given by Eq. ( 2.12| ) 
appears naturally in the theory. The basic averaging equation is Eq. ( p. 11 ), which enables the 
computation of averages of functions that depend on the sum of the squares of the interaction 
potentials in terms of v. This is equivalent to the distribution of Eq. (2.23) for the quantity = 



Sfc where Vk is the dipole-dipole interaction potential defined by Eq. (2.2). The initial rate of 
conversion from ss' to pp' is given by Eq. ( 2.37 ) or Eq. ( 2.85| ) to be ^/ttv/2, and is independent of the 
detuning A. As we will see later, this initial rate remains the same even when the sp-ps interaction 
is switched on. 



Chapter 3 



Localization in the Simple Model 



3.1 Introduction 



In this chapter we continue with our discussion of the simple model of Chapter ^ by examining 
the issue of localization. We want to determine the effective range, or localization length, of the 
disturbance caused by the presence of an s' atom in the frozen gas of s atoms. Since the 
interaction does not have a length scale, it follows that the localization length must be of the form 
(n/N)^^^ g {vt, A/v), where N/il is the density of s atoms, t is the time, A is the detuning from 
resonance, and v is defined in Eq. or more generally in Eq. (2.73). The task, then, is to give 

a precise definition of the localization length and to compute the function g (vt, A/w). 



First, in Section 3.2 we compute a distribution function (5(A, r, i)), which describes the con- 
tribution to the ss' —>■ pp' process due to s atoms at position r when the s' atom is at the origin. 
We begin by computing the Laplace transform of this function in Section [3.2.1 , then we invert 



this Laplace transform in Section |3.2.2[ In Section |3.2.3| we perform a check on the results of Sec- 
tion 



3.2.1 and 3.2.2 by using them to rederive the result for the signal that was found in Chapter 



In Section 3.2.4 we present plots of the distribution (5(A,r,t)). 



In Section 3.3 we compute the moments of the distribution (5(A,r, i)). We again start by 



computing the Laplace transform of the moments in Section [3.3.1 , then we invert this Laplace 
transform in Section 3.3.2[ . The result of Section 3.3.2 is rearranged in Section 3.3.3 to obtain a more 
convenient form, and in this section we also define localization lengths in terms of the moments in 



the usual way. The end result for the moments is checked in Section 3.3.4 by seeing that the zeroth 
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moment agrees with the resuh for the signal as a function of time that was found in Chapter ^ In 



Section 3.3.5 we present plots of the results. 



The summary and conclusions are given in Section 3.4. It is possible to go to this section directly 



and refer back to final formulas and the plots of Section 3.3.5, skipping the intervening mathematics. 



3.2 The function {S (A, r, t)) 

In order to study localization in the simple model of Chapter |[ we consider the quantity 



N 



(5(A,r,t))^(^|c,(A,t)|^5(r-r,)). (3.1) 



As is illustrated explicitly below, the computation of the average is very similar to the computation 
of the average signal, (S* (A, a)), of Chapter p|. 



3.2.1 Computation of the Laplace transform of (5(A,r,t)) 



From Eq. (|2j) of this thesis or Eq. (3) of Ref. ||l§ we have that 



where — J2k^k- Therefore we define 

Consider for now just the fc = 1 term in the sum. For convenience, we will call this term A. We 
have 

o^^"" J a2 + A2+4F2(ri) + 4^;i2'>4 

where I have written Vi as V (ri) to be explicit. The integral over ri is trivial, and we find 

A = -V^ (r) [ d^r2 ■ ■ ■ d^rN ^ • (3.5) 

a ^ J a2 + A2+4V2(r)+4E/l2"^^ 

We observe now that we can write the denominator as an integral over an exponential and use the 
same averaging technique as in Chapter |2[ The fact that the sum over I starts from I — 2 does not 
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affect the validity of Eq. (2.11), as long as N is large. Thus we have 



dA / (fr2 ■ ■ ■ d rjv exp ■! -A 



N 



(r) +4^V;^ 



2 1 f r -1 

-V^ (r) - J dAexp |-A [a^ + + AV^ (r)] - 2wVa| 



(3.6) 



It is clear that it did not matter that we chose the fc = 1 term from the sum in Eq. (3.3). We 
would have gotten the result of Eq. ( [3.6[ ) no matter which term we picked. It follows that all the 
terms are equal, and so 



S (A, r, a) ) = - — V^ (r) J dA exp |-A [a^ + A^ + AV^ (r)] - 2wVa| 



(3.7) 



The remaining integral over A is familiar from Chapter ^ and can be looked up in tables (see, for 
example, Rcf. p^) or done by Mathematica. The result is 



5(A,r,a) 



N 



21/2 (r) 



2^05'^^ (r) 



exp 



a [a2 + A2 + 41/2 (j.)] 
i>2 

_a2 + A2 + 41/2 (r) 



3/2 

erfc 



v/a2 + A2 + 41/2 (r) 



(3.8) 



3.2.2 Inverting the Laplace transform 



Now that we have (S (A, r, a) ), the next step is to invert the Laplace transform. Because Eq. (3.S) 



is of a form that is similar to Eq. ( 2.27 ), we can again use Eqs. ( 2.31 ), ( 2.33 ), and ( 2.34 ) to achieve 
this. 



Using Eq. ( 2.31 ), we have 



5(A,r,a) 



N 2T/2 (r) 

a [a2 + A2 + 4^2 (j.)] 

N 2v^V^ (r) 



r)f/M^nr(f + 1) 



a [a2 + A2 + 41/2 (r)]=^/' t'o ^ (t + l) 
N 2T/2 (j.) 

a [a2 + A2 + 41/2 (j.)] ~ 

?7 



^^0,2 + A2 + 41/2 (r) 



^2.V^F2(r)^ ^ 



-vY 



r (t + 1) a [a2 + A2 + 4y2 (r)]"^ 



(3.9) 
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Next, using the fact (from Eq. (2.66) of Ref. pM) that the inverse Laplace transform of Eq. (2.32) 



is Eq. (2.34), we see that 



(5(A,r,i)) 



N 

n 

N 

n 



2V'^(r)^iF,|l;|2;-i[A2 + 4y^(r)]t^ 
2vV^V' (r) y / X 

%eor(f + i)r(n + 4) 

n + 3 n + 4 n + 5 1 [^2 + 4^2 (^^^ ,2 



^ 2 ' 2 ' 2 ' 4 

N 

n 

^2V^vV' (r) t' f -^-^ 

^ ^„ r (f + 1) r (n + 4) 



;i/^(r)tSF2|l;|2;-i[A^+4V^(r)]t^| 



1F2 



ri=0 

n + 3 n + 4 + 5 1 



2 ' 2 



[A2 + 4^2 (r)] t2 



3.2.3 A check on the result for (S ( A, r, a 



It is clear from Eq. (3.1) that 



(3.10) 



(5(A,r,t)) = (^(A,t)), 



(3.11) 



and since taking the Laplace transform with respect to time and integrating over the variable r 
clearly commute, it must also be the case that 



d'r 5(A,r,a)) = (S{A,a) 



(3.12) 



For the sake of completeness, we now check that this is true. We will also find this exercise useful 



in the next section. From Eq. (3.7) we have 



rf-^r 5(A,r,«) = 



,3 AT 21/2 (r^ 



n a 

dA cxp 



a JO 

1\ _d_ 
4 j OA 



dA exp 

N 



dAexp|-A[a2 + A2+4F^(r)] -2u\/a| 
1 N f 

-A {a^ + A^) - 2i;%/Aj — / d^r (r) exp [-AAV^ (r)] 
A {a^ + A2) - 2t;VAj x 
d^r exp {~4AV^ (r)) . (3.13) 
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We now recall from Eqs. (2.10) and (2.12), or more generally from Section 2.5 of this thesis, that 



N 



It follows from this that 



and so 



N 

n 



J (fr exp [-4Ay2 (r)] = iV _ 2wVA, 



(3.14) 



(3.15) 



d^r (5(A, r, a)) = — / dA exp 



2a 



dK 
2a Jo 



exp 



-A (a^ + A^) - 2wVA 
-A ia^ + A^) - 2v^/^ 



(3.16) 



The integral over A can again be found in tables or done by Mathematica. The result is 



dV(5(A,r,a) 



'TT V 



2 aVa^ + A2 
S{A,a)), 



exp 



erfc 



Va^ + A^ 



(3.17) 



where we have used Eq. ( |2.27D . 

A similar check for {S (A, r, t)) is carried out in Section ^.3.4 



3.2.4 Plots 

In Fig. [3.1| we see plots of the dimensionless quantity 

r2 (5(A,r,i-oo)) 



/ \ 2/3 r> 



\flfl'J N 



(3.18) 



1 /3 

as a function of {v/iifi') ' r, for A = 0, v, and 2v. The quantity (5 (A, r, i = oo)) is most easily 



computed by multiplying Eq. (3.8) by a and then taking the limit a — )■ 0. We see that the peak of 
the distribution shifts to the left with increasing A. 
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Figure 3.1: Plots of the quantity (5 (A, r, t = oo)) as a function of r for A equal to 0, v, and 2v. 
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3.3 The moments of {S (A, r, t)) 



In order to compute the moments of {S (A, r, i)), we first define 



(3.19) 



so that A4p{A,a) is the pth moment of (^S {A,r, a)j . We will see that we can then determine 
Mp (A,t) by applying the inverse Laplace transform methods we have already refined. 



3.3.1 Computation of Alp (A, a) 



Much of this computation is very similar to that done in Section 3.2.3. Using Eq. (|3.7D, we have 



MpiA,a) 



^3 „N2V^r) 



n a 

dA exp 



dA exp [-A {a^ + A^ + (r)) - 2vVA 



dA exp 



a JO 
2 

a .10 

l_d \ (N 
AdAJ 1 n 



-A {a^ + A^) - 2uVA 
-A {a^ + A^) - 2vVA 
J d^rrP exp [-AAV^ (r)] 



N 



d^rrPy^ (r) exp [~AAV^ (r)] 



(3.20) 



We now note that 



J d^rrPexp [-iAV'^ (r)] = J J '^^^^^ {l - exp [-iAV'^ (r)] } 



(3.21) 



Although the first term diverges, it will disappear when we take the derivative over A because it 
does not depend on that variable. 

We now need to compute the remaining integral over r. I will compute it here for the case 
where V (r) — fiji' jr^. If there is an angular dependence in the potential, then the computation is 
somewhat more complicated because the extra factors of r in the integrand will bring down extra 



angular factors when we apply the methods of Section 2.5. This is one case where we cannot simply 



rescalc v to obtain the correct result when an angular dependence is introduced into the potential. 
With this in mind, we have 



j d^rrP {I- exp [-4AF2 (r)] } = At: j dr rP+^ J 1 



exp 



-4A 
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47r 

T 



4A(mm')^ 



dz zP/^ 



1 — exp 



, (3.22) 



where = rV4A(/i/i') ■ The remaining integral over z can be found in tables or done by Mathe- 
matica. The result is 



dVrP {1 - exp (r)] } = -\r (-^) f [4A i^.f^'y 



(3.23) 



for p < 3. This should sufhce, since we will be primarily interested in p = 1 and p = 2. 
Thus we have (for p < 3) 



MpiA,a) 



1 f°° 

— / dA exp 
2a io 



-A(a^ + A^) -2vVA 



1, 



dA 1 2 
1 

^2^ 



p + 3\ 47r 

P + 3 



4A(mm')^ 



X 

p + 3 



4(a^m')' 



— — / dAA 6 exp 



-A {a^ + A2) 



2vy/A 



p-3 

dAA 6 exp 



'A (a^ + A^) - 2wVA 



(3.24) 



where in the last step we have used the fact that F (n + 1) = nT (n) and Eq. (2.12). For p > —3, we 
have from Mathematica that 



p-3 

dAA 6 exp 



'AA^ - 2v^/A 



2v ^ + ^ fp + ^. 3. 



. P + 6 

A— 



6 '2'A2 



(3.25) 



where as in Chapter ||, A^ — a'^ + A^. It follows that for \p\ < 3, we have 



a2 + A2 
2v 



6 



6 '2'a2 + A2 



(a2 + A2) 



6 ; ' ' ^ 6 '2'a2 + A2y' 



(3.26) 
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3.3.2 Inverting the Laplace transform 



We know from Eq. ( p.35D that 



iFi (a, b, z) 



r (j + a) 



£(6) 

r(a)^r(j + 6)r(j + i)- 



(3.27) 



so 



1 



Q,2 + A2 

2?; 



P + 6 

(a2 + A2) — 



p;r(j + i)ro- + i)(^2 + A2)^+* 



^v^E 



,t^r(j + |)r(j + i) («2 + A2)^+ 



P+6 
6 



E 



,2i 



E 



J 

2] + l 



^r(j + i)r(j + i) («2 + A2)^+ 



P + 3 
6 



^r(j + |)r(j + i) («2 + A2) 



(3.28) 



where we have used the facts that F (1/2) = -^7? and F (3/2) = Vtt/^- 
Now we note that 



^ F0- + H±3) ,2, ^ r(, + i^) 

r (j + 1) F {j + 1) + A2).+4^ F (, + 1) r (j + 1) + A2)^-+4^ ^ ■ ^ 

can be rewritten as 



^or(¥ + ^)r(¥ + i) V V^^^+A^; V«^ + A2 
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3.3.3 A more convenient form for Aip{A,t) 



As we mentioned in Section 2.4.2 , it is sometimes useful to have a result as a series in powers of the 
detuning squared. Now we exploit the method used there to obtain such a form for A4p (A, t). For 
convenience, we define 



r(^ 



^r(i±i)r(i±2)r(, + 2 + |) 



1F2 



i + 1 p, j + 2 p i + 'i 

2 6' 2 6' 2 



(3.36) 



so that 



M, (A, t) = 2^r f i-^) (MM'f vt'^B 



We now recall from Eq. ( 2.35 ) that 



(3.37) 



(a; 6, c;z) 



r (6) r (c) ^ 



r (fc + a) 



r(a) ^r(fc + 6)r(fc + c)r(fc + i) 



(3.38) 



and hence 



E 



^r(i±i)r(i±2)r(, + 2 + |) 
r(fc + ^ + f) 



(-wi)-' X 



^or(fc + ^ + i)r(fc + ^ + f)r(fc + i) 

Recahing Eq. ( |23o| ), we see that 



(3.39) 



(3.40) 



and so 



00 00 



B = 



j+3 E/ 



£+3 ^ ^ p / j-^i ^ p 



r(fc 



i±i 
2 



r (fc + + f ) r (A; + ^ + f ) r (fc + 1) 



--A^t^ 
4 



r(fc + i) 
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^^„r(4i)r(42)r(i±^ + ^. + |)r(43 + fc + |) 1, 2 



Now we define 



K?odd 

' ",.otoddr(^)r(^)r(^ + fc + |)r(^ + fc + f) V 2 

and similarly, 

r(4i + fc + f) f vt 



K7CVCn „ \ 

' ",.oreve„r(^)r(^)r(i±^ + fc + f)r(^ + fc + |) ^2 



j—0,j even 

SO that 



B = #|(-^)'niTi)(''r + «r") 



We then have (again using Eq. ( 2.35| )) 



2 ^^^ro + i)r0- + i)r(j + fc + | + |)r(j + fc + | + 2) 

r(fc + | + i) 
2r(|)r(fc + f + |)r(fc + | + 2) 

6 '2' 6 2' 6 ' 4 
1 



V^(fc + | + l)r(fc + | + |) 



n 3 p 3 n w2^2 



Performing the analogous calculation for SS™", we see that 



.even _ T (j + + f + i • 



CX2 

6' 



,t^r(j + i)r(, + i)rO- + fc + | + i)rO- + fc + | + |) 



r(fc + f + ^) 



r(i)r(fc + f + i)r(fc + f + |) 

/, n 1 1 , p , p 3 v'^t^ 

1F3 fc + - + -; -,/c + - + 1,A: + - + -: 

6 2' 2' 6 ' 6 2' 4 

1 1 



V^(fc + | + i)r(fc + | + i) 

/, p 1 1 , p , p 3 -y^t^ 

^ H 6 2' 2' 6 ' 6 2' 4 
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Using Eq. (3.44), we arrive at 



B = 



^ OO 



2^ 



fc=0 



r(fc + i) 



(fc + | + i)r(fc + | + |) 



F. k 



1-^3 



(fc + | + i)r(fc + | + i) 



3 , p 3 , p u^i 

l:-,/j+- + -,fc + - + 2: 

'2' 6 2' 6 ' 4 



2+2 \ n 



(3.47) 



Recalling Eq. (3.37), we obtain finally 



Mp{A,t) 



2 3 r I — I i^ifi ) Vt 3 X 



i±3 $Z 



24^ 



fc=0 



4 

1 1 



r(fc + i) 



^^^1^ + 1+2^2'^ + ! + ^'^ ■ 6 ■ 2' 4 



(fc + | + i)r(fc + | + i) 

p 3 v'^t'^' 



(fc + £ + i)r(fc + | + |) 



1^3 A: 



6 



4 I 6 



E 

fc=0 



A2t2 



1 3 , p 3 , p 
1; -. fc + - + -, A: + - 

'2' 6 2' 6 

— (wi) 3 X 

V J 



1 



2+2 \ 1 



r(fc + i) 



(fc + |+i)r(fc + | + i) 

2+2 



1^ I , P 11. P 1, P 3 V t 

1F3 A; + - + -; -, fc + - + 1, fc + - + -; 

6 2' 2' 6 6 2' 4 



(fc+| + l)r(fc 



iF3(^ + | + l;-,fc + f + -,^ + f + 2; — 



(3.48) 



It is important to recall that this result is valid for \p\ < 3. 

The result of Eq. ( [3.48| ) shows that the pth moment A4p is equal to the quantity {fifi' /v)^^'^ 
multiplied by a dimensionless function. Since pp' /v is proportional to Q/N, as can be seen from 



Eq. (2.12), it follows that Mp has the dimensions of length to the pth power as it should. 

To extract a localization length Ip from the moment Mp, we must first normalize the distribution 
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S of Eq. (3.1C). The localization length corresponding to the pth moment is thus given by 



lp{A,t) 



Mp{A,t) 
L-A^o(A,t) 



i/p 



(3.49) 



For each value of p, this expression gives a somewhat different measurement of localization. The 
quantity Ip is proportional to (/^/i'/w)^^^, which according to Eq. ( 2.12| ) is equal to (3/47r'^/^)"'^^'^ (fl/N)^^'^, 

1 /3 

or 0.5126 i^/N) ' . Therefore we see that Ip has units of length for all values of p, as one would 
expect. 



3.3.4 A check on the result for Mp{A,t) 



We see from Eq. (3.11) that Mq t) = {S (A, t)), and this provides one way to check Eq. (3 



From Eq. ( 3.48 ) we have 

1 /I \ °° 



1 



(fc + i)r(fc + i)^^n' + ^^i' + ^'' + l^' 



vt 



{k + l)T{k 



3 3 v'^t'^ 
-1F3 ( fc + l;-,fc + -,fc + 2; — 



fc=0 



vt 



11 3 v'^t^ 

.2-1^3 ( fc + -;-,fc + l,fc + 



2' 2' 



F (A: + f ) F (fc + 2; 



3 3 
iF3[k + l;-,k+ -,k + 2: 



2' 4 



which agrees perfectly with Eq. (|2.4l| ). 



(3.50) 



3.3.5 Plots 



Figs. 3.2 - 3.4 give plots of the dimcnsionlcss quantities 

p/3 



Mp{A,t) 



(3.51) 



as a function of vt, for p = —1, —1, 1, and 2 and for A — 0, v, and 2v, respectively. 
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In Figs. 3.5 - 3.7 we see plots of the dimensionless quantities 

1/3 



(^) '^<^'"> 



(3.52) 



as a function of vt, for A = 0, w, and 2v, respectively. We see that all the plots initially rise 

1/3 

rapidly like {vt) , and then oscillate around a value of order one that depends rather weakly on 
the detuning A. 



3.4 Summary 

We study localization in this system because we need to have some idea of the relative importance 
of close pairs of atoms, as well as some idea of the effective range of the interaction. Naively, it 
seems that the close pairs are dominant and that the range of the 1 /r^ potential is infinite because 
diverges at r = and yet the integral J d?r/r^ over all space diverges at both limits. 
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Figure 3.2: Plots of the moments Aip {0,t) for p equal to —2, —1, 1, and 2. 
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1.5 




Figure 3.3: Plots of the moments Mp {v, t) for p equal to —2, —1, 1, and 2. 
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Figure 3.4: Plots of the moments Mp {2v, t) for p equal to —2, —1, 1, and 2. 
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Figure 3.5: Plots of the localization lengths Ip {0,t) for p equal to —2, —1, 1, and 2. 
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Figure 3.6: Plots of the localization lengths Ip {v,t) for p equal to —2, —1, 1, and 2. 
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Figure 3.7: Plots of the localization lengths Ip {2v,t) for p equal to —2, —1, 1, and 2. 



Chapter 4 



Sparse ss^ pp^ Processes with U 
Interaction 



4.1 Introduction 



In this section we complicate the simple model of Chapter g by introducing the sp — s- ps process 
through the addition of an interaction potential U. The sp ~^ ps process can be viewed as an exciton 
propagating through the medium |||, |l^, |3^ , as was discussed in Ref. . 



In Section 4.2 we discuss the Cayley tree approximation, starting with its justification in Sec- 



tion 



4.2.1. In Sections 4.2.2 - 4.2.4 we solve for the Laplace transform of the probability amplitude 



in the case where A and V are both zero, the case where V is nonzero but A is still zero, and the 
case where both A and V are nonzero, respectively. 



We consider a Hubbard band approximation in Section 4.3. In Section 4.3.1 we again solve for 



the Laplace transform of the probability amplitude in the case where A and V are both zero and in 



Section 4.3.2 we do the same for the case where V is nonzero but A is still zero. However, we note 



in this case that the latter solution is not consistent. In the case of the Hubbard approximation, we 
are also able to compute the square of the probability amplitude in the V — A — case, as is done 



in Section 4.3.3 



In Section 4.4 we consider a coherent potential approximation (CPA). In this approximation we 



are able to compute the probability amplitude in all cases. We will see in Section 4.6 that the CPA 
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results agree best with the resuhs of the numerical simulations discussed in Chapter 0. 



In Section 4.5 we compare the results of these three approximations. We see that they are in 
some ways rather similar, but that the differences are sufficiently large that it is easy to tell them 
apart, particularly when the sp — > ps process is weak. This means that one approximation will 
fit the exact results better than the others, and in the case we are studying it is the CPA. Thus, 



in Section 4.6 we present plots of the CPA results in comparison with the results of the numerical 



simulations discussed in Chapter 0. 



Finally, in Section 4.7 we summarize these results and conclude. 



4.2 The Cayley tree approximation 

4.2.1 Justification of the Cayley tree approximation 

We start with the equations 



N 

ido = Aao + ^VfcCfc, (4.1a) 

k=l 

ick = Vkao+^UkiCi. (4-lb) 

l^k 



These are the same as Eqs. ( 2.1a ) and ( ^.Ib ), except that we now allow for the sp ps process 



by including the Ukici term in the equation for Cfe. For the purposes of compressing indices on 
sums and products, it is convenient to define Uki to be zero when k — I. With this definition, the 
previous equation becomes 

N 

iao = Aao + ^VfcCfe, (4.2a) 
fe=i 

N 

iCk = Vkao + ^UkiCi. (4-2b) 

1=1 

Unfortunately the potential U makes the set of equations insoluble as they stand, so we cannot 
obtain an exact solution as we did in Chapter ^. However, the equations can be solved numerically, 
and one can also make quite a lot of progress by making approximations and proceeding analytically. 
In developing approximations, it is useful to cast the problem in the standard Green's function 



language ||l^ , as discussed in Section 3 of Ref . [gT[ . The Green's function of Eq. (|4.2b[ ) satisfies the 
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equation 

N 

^Gin — Sin + ^ UimGynn, (4.3) 
m— 1 

and then 

ao (A, ^) . ■ (4.4) 

Actually, — is the 00 element of the Green's function for the entire set of + 1 equations. 



We can derive a convenient form for the on-site Green's function G„„. From Eq. (4.3) we have, 
for I ~ n, 

^Gnn = 1 + UnmGmm (4.5) 
m 

and if we write separately the term containing G„„, we obtain 

wG;„ = UlnGnn + ^ UlmGmn, (4.6) 



for I n. We note here that Eqs. (4^) and ( ^ ) are the analogs of Eqs. ( 4.2a ) and ( 4.2b| ), 
respectively. 

Now we imagine removing the nth atom. The new Green's function for the remaining — 1 
atoms, which we will call Gip_[„], obeys the equation 

"^Gjpjn] = ^ip + ^ UlmGjnp^[ri], (4.7) 

where / and p are never equal to n. If we multiply both sides of this equation by UpnGnn and sum 
over p, we obtain 

] Glp,\n]UpnGnn — UlnGnn + ^ ] t^im G,„p_ [„] ?7p„G„„ . (4.8) 



Comparing this equation with Eq. (|4.6| ), we see that for m ^ n 

Gni7i — ^ Gjyip,\n]UpnGnn- (4.9) 
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If we insert this result into Eq. (4.5), then we have finaUy that 



and hence 

It is important to keep in mind that these Green's function equations are exact. 

The Cayley approximation ^ keeps only the diagonal {I — m) terms in the sum over 

I and m in Eqs. (14) and ( 4.1l| ). It further assumes that each G;; [„] is independently distributed, 



and that there are sufficiently many sites that Gu,[n] and Gu have the same distribution. 

We choose to set lj — ia and work in Laplace space instead of Fourier space because this makes 
many of the integrals that follow more obviously convergent. Up to a point, the same analysis can 
be carried out with the Fourier analogs, but Laplace transforms are needed to solve the Cayley 
equations exactly in Section pi. 2. 2 . We also work with the inverses of these Green's functions, setting 



ao (ia) = l//o (o?) and Gkk {ioi) — i/ fk (a), so in the Cayley approximation we have from Eqs. (4.7) 



and (4.8), respectively 



T/2 

/o = a + zA + ^-p, (4.12a) 

k=i 

Ik = « + Et^- (4-12b) 

1=1 



4.2.2 Cayley tree approximation with the U process only 



First we consider just the U process and solve Eq. (4.12b), which is to say we look at the sp — > ps 
process by itself. Note that fk is always real and positive. The distribution of fk is given by 



P{fk)= / {5\fk 



N TJ2\\ N 

^1) n ^(/™)^/" 



(4.13) 
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where for a uniform gas of volume fi, 



N 



(4.14) 



Writing the delta function as a Fourier transform, we see that 



P{h) 



igifk~a) 



N 

J n 



exp ~iq 



\ " ^ 



N 



( U'^ 

I 11 d'riexpl-iq^ 
■' i^k,i=i ^ ■' 



Using the fact that (J]™ Pm) (Hi Ei) = Dm PmEm, we can write 



(4.15) 



Pifk) = 



r p. e^«(/.-) f[ [ P (/„,) 1 / d3 r„ exp 



-iq 



Mifk~a) 



— OO 
OO 



dq 
27r 



If f Ip 



fi 



TP 



dfPif)^ J d\ |l-exp 



-iq- 



f 



(4.16) 



where in the last step we have simply changed variables from ri to r = ri — r^. Using the averaging 
method detailed in Section 2.3, we see that]] 



dq 
2^ 



Pif)^ I ^ e'-J^-^-") exp 



u\JiqQ (a) 



(4.17) 



where 



and 



Q(a)=^°°-^P(/) 



47r3/2 



(4.18) 



(4.19) 



^Of course, we must use the more elaborate averaging method of Section 2.5 if the potential V is more complicated 
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if U is of the simple p? /r^ form with no angular dependence, and is the equivalent of v in Eq. (2.73) 
otherwise. For A and B real, we have that 



dq ^qA^^qB _J_ ^ (-—\ e (A) 



(4.20) 



where this integral is evaluated by the same procedure used to evaluate the integral of Eq. (2.15).^ 
With this result, we have 



and so 



We also have 



P{f) = 



'{a) 



uQ (a) 



2^ i^j _ afn 



df 



exp 



u^Q^ (a) 
4(/-«) 



'(/-«). 



uQ (a) 



— a 



3/2 



exp 



u^Q^ {a) 
4 (/-a) 



: exp 



u^Q^ (a) 



4a 



erfc 



uQ (a) 



2^ 



(4.21) 



(4.22) 





/•OO 




Jo 




1 




a 



/tF uQ (a) 
2 a3/2 



exp 



u^Q^ (a) 



4a 



erfc 



uQ (a) 



(4.23) 



or 



1 

a 



(4.24) 

Thus we can solve numerically the transcendental relation for Q (a) given in Eq. (4.22) and plug 



it into Eq. ( 4.24 ) to find (1//). It is important to note that Eq. ( 4.22 ) defines Q as an analytic 
function of a in the half-plane where Rca > 0. Thus (1//) can be analytically continued, in the 
variable to — ia, from the positive imaginary axis to the entire upper half-plane. 



^In case the reader is curious, the thcta function appears here but is omitted in Eq. (2.23) because is an 
intrinsically positive quantity. 
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4.2.3 Cayley tree approximation with U and V processes in the absence 
of A 

Now we allow the V process to be present, but we still ignore the detuning, A. We start with the 
equations 



N 



k=i -"^ 



fk = a+ 



l=l,l^k 



^kl 
fl 



(4.25a) 
(4.25b) 



Repeating the analysis we did before, we have 



Poifo)= J /6(fo-c.-f2j^))l[Pifi)dfi 



k=i ■' / / 1=1 



(4.26) 



which leads to the form 



^0 (/o) 



dq 
2^ 



MUa-oi) 



-iq- 



(r) 



/ 



N 



(4.27) 



At this point it is clear from the work we did to obtain Eq. ( 1.21 ) from ( 4.16|) that this equation 
leads to 

v^Q^ [a) 



4(/o-«) 



(4.28) 



where Q (a) is still given by Eq. ( [4.22|) , and depends on u. Therefore we obtain the result 



rfpoifo) 

Jo Jo 

1 ^/T:vQ{a) 
a ^ ~2 ^^372" 



■ exp 



4a 



erfc 



vQ (a) 



(4.29) 



which again can be analytically continued. 
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4.2.4 Cayley tree approximation with U and V processes in the presence 
of A 



If we further consider the effect of a detuning, A, then we have the fuU set of Eqs. ( 4.12a ) and 
( 1.12h| ). The equation for fk is unchanged from the previous case, but because the equation for /o 
is now complex we have to consider a distribution in the two variables — Re /q and /g = Im /q . 
We have 

^o(/o^/oO= / (5[^f^~a-Y,^y{f,-A)^\{PU,)df,. (4.30) 

Clearly the delta function involving /q does not participate in the averaging or in the integrations 
over Therefore we can treat the part involving just as we did before, and the /q part will 
just be carried along. Hence we find 



vQ (a) 



20F(^H 



\3/2 



exp 



v'^Q'' (g) 
'4(/o^-a) 



{l^-a)6{f,-A). (4.31) 



It follows, then, that 



^■^0 fR , . fi Pa (/(f ) /o) 

) JQ + ''JO 



1 



vQ (a) 



a + iA 



{a + iA) 



3/2 



exp 



(a) 



4(a + iA) 



erfc 



vQ (a) 
+ zA, 



(4.32) 



4.3 The Hubbard approximation 

It is also possible to make approximations other than the one we have chosen. For example, instead 
of working with Eq. ( 4.12bD , one can instead choose to work with 



AT 



1=1 



which gives explicitly 



(4.33) 



(4.34) 



This approximation corresponds to a Hubbard band. 
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4.3.1 Computation of (1//) 



It is possible to average Eq. ( 4.34 ) , although it requires a little more manipulation than the averaging 
methods we have applied previously because of the presence of the square root. Wc first note that 



1 

fk 



a + \Ja^ + 4ZY2 



(4.35) 



where 



(4.36) 



Interestingly enough, the Laplace transform in Eq. (4.35) has a known inverse. We see from 
Eq. (2.1.35) of Rcf. that the inverse Laplace transform of Eq. (|4.35|) is 



-.Km)- 



(4.37) 



We will use Eq. (4.37) in Section 4.3.3 to compute ^|ao (A — 0, i)| y in the case where only the \J 
process is present. 
We now note that 

y-x^^ , ^ , (4.38) 



2 Jo ^/a^+Xy 



so we can rewrite Eq. (4.35) as 



1 



fk Jo Va^ + 4AZ^2 
We introduce yet another integral by observing that 



if Re (CA) > 0, so 



1 2 r°° 

— = ^ / dp dXexp (a^ + 4AZY2)1 
fk Jo Jo 



From our previous experience with averaging over U, we know that 



j'^ = -^J^ d/3 dX 6x^(^-/3"^ + 2upVX 



(4.39) 



(4.40) 



(4.41) 



(4.42) 
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Here we write (1//) in place of {i/fk) because the averages of all the fk are equal. The integral 
over f3 is familiar to us, since it is the type of integral that shows up repeatedly in this averaging 
method. We have 







- / 






Jo 



2a ' 



OA exp — ^ 
V " 



erfc 



(4.43) 



The remaining integral over A can be done by parts, with the result that 



U\/TT U 



exp 



erfc 



- 1 



(4.44) 



We can also compute (1//) by first computing the distribution Pjj (/fe)- It will turn out that 
knowing Ph (fk) will aid greatly in the computation of l//o for this approximation, as we will see in 
Section 4.3.2| . Because in this case fk is explicitly given in terms ofU^ through Eq. (4.34), we have 



a 1 



(4.45) 



This average is easy to compute using the analog of Eq. (2.23) for the potential U. We have 



/•OO / -1 

1 2/ - a 



[/(/-«)f/' 
where we have used the fact that 



exp 



4/ (/ - a) 



1 u 
9{f-a) 



exp 



u 

lz72 



(4.46) 



a 1 

2 ~ 2 



Va^ + j = (2/ ~a)S[U'~fif- a)] 



(4.47) 



The theta function is present in Eq. (4.46) because and / are positive definite quantities, and 
hence S [U'^ ^ f if ^ «)] cannot be nonzero unless / — a > 0. 



Using the distribution of Eq. ( [4.46D we can compute 

df 



/ 
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df I 2/ - g 
1 r df 2f~a 



exp 



2V^7a / [/(/-a)] 



3/2 



exp 



4/ (/-a). 
"4/ (/-«) 



'(/-a) 



(4.48) 



This integral is very difficult to compute analytically in its present form, even by a computer package 
such as Mathematica. Therefore we change variables to % = /(/ — a). Then dx = (2/ — a) df, and 



because / must be positive we have f — a + ^/a^ + Ax. It follows that 



dx 



f I 2^ Jo x^''^ a + yja'^ + 4x 



exp 



4x 



(4.49) 



This integral can now be done by Mathematica, and it gives exactly the result of Eq. (4.44) 



4.3.2 Computation of (l//o) 



Going back to Eq. (4.12a), we see that 



/o 



(4.50) 



Since this is the same equation we considered for the Cayley approximation in Section 4.2.4 we see 



that the result for Pq.h (/o^, /o) is very similar to Eq. ( 4.31 ). Specifically, we have 



-Po.H (/(f , /o ) 



1 vQYi (a) 
(/«-«)'/' 



exp 



v-'Ql (a) 



where 



4 {!i 



df 



{fo^~a)5{f^~A) 



It then follows from Eq. (4.32) that 



(4.51) 



(4.52) 



1 



a + ilS. 



/tF wQh (a) 
2 (a + ?A) 



3/2 



exp 



(a) 



4(a + iA) 



erfc 



«Qh (g) 
2Va + zA 



(4.53) 



One might think that because Pjj (/) is known explicitly in this case, instead of involving tran- 
scendental quantities as Eq. ( 4.21[ ) does, it should be possible to compute the function Qh ipt) 



explicitly as well. Unfortunately the integral of Eq. (4.52) cannot be done analytically. If we make 
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the same change of variable that we did to obtain Eq. ( 4.49| ), however, then Eq. (4.52) becomes 



Oh (a) 



(4.54) 



and this integral can be done numerically for a given value of a. 



Thus Qh (q;) can be computed numerically, and the result substituted into Eq. (4.53) to obtain 
(l//o). 



Although we have shown how to compute (l//o) for the Hubbard model using Eqs. (4.12a) and 



(4.34), we note that these equations are not consistent. Consistency demands that (l//o) reduces 



to (1//) when u ~ v and A = 0, and this is clearly not the case with Eqs. (4.53) and (4.33). In 
fact, one can see this inconsistency immediately from Eqs. (4.12a) and (4.34). If we set if' = 
and A = in these equations, then we have 



/o 
/fe 



T/2 



N 

/fc ' 

fk 



(4.55a) 
(4.55b) 



We see that even in this limit, /o and fk are still treated differently in this approximation. Therefore 
the Hubbard approximation is inconsistent, and to be safe one should disregard the result for (l//o) 



and only consider the result for (1//) of Section 4.3.1 as valid 



4.3.3 Computation of (^\ao{A,t)f) with the U process only at resonance 



We saw in Eq. ( 4.37 ) that at resonance 



ao(0,i) = —Ji{2Ut), 



(4.56) 



and so 



|ao (0,t)r 



Jl {2Ut) . 



According to Eq. (10.1.1) of Ref pW, this function has the inverse Laplace transform 



- -S'(0,a) 
a 



— - + 2 - 2 cos 6* 

4W2 2U 



(4.57) 



Chapter 4: Sparse ss' — + pp' Processes with U Interaction 



67 



de{l + cos0) ^/a^ + 8U^{l-cose)-a 



(4.58) 



Now we again apply Eq. (|4.38| ), with the result that 

■ / d9 {1 + cos 9) 4 {I- cose) [ 
Jo Jo 



1 1 

1 2 r 



dX 



y/a^ + 8\W (1 - 



d\ 



de SIX? 9 , , 

a TT ./o Jq + 8AZ^2 _ ^os ( 



At this point, we apply Eq. ( [4.4C| ) to find that 



(4.59) 



1 4 

S{0,a) = - 



a TT 



3/2 







[ d9 sin^ 9 


fdxj 


'o 


Jo Jo 



-/3 a^+8AW^(l-cos« 



(4.60) 



This expression can now be averaged. We find that 



1 4 

S{Q,a)) = - 



a TT' 



3/2 



/•TT 




/ d9 sin^ e 


fdx 1 


lo 


Jo Jo 



(4.61) 



The integration over /3 can be performed at this point, with the result that 

4 



S{0,a) 



1 

a 

exp 



d9 sin^ 



3/2 

2 (1 — C0S( 



TT-' - JO 





erfc 



dA -^^ X 

2 a 



a 

/ d9 sin^ 9 dX exp 

a TT a Jo Jo 

erfc 



V2A(1 - cos( 



2Am2 

— ^ (1 - cos( 



-V2A(T^^^cosl 



(4.62) 



The A integral can be done by parts. The result is 
'S{0,a) 



a 



I 21 r 9 

/ d9 sin L , 

ct TT a ,/o [2m^(1 — cos( 



■ exp 



2m2 

^(l-cos6') 



erfc 



u 



\/2(l-cos( 



a 



^-^4 / 



sm 



1 — cos 



2^2 (1 — cos 
- exp 



2^2 

(1 - cos( 



TT UVl— cos 

u 



f d9 (1 + COS0)- I - 

TTU^ Jo 



2x3/2 ^ 



TT / U 



sm' 



erfc 

2 /) 



-a/2(1-cos( 



d9- , 

vl — cos( 



(4.63) 
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The last two terms are standard integrals, and we have that 
' S{Q,a) 



1 la 

/ ^6* (1 + cos6')exp 







a TT u 

a 16 1 
~ ~ 37r3/2 - 



^(l-cos0) 



erfc 



-v/2(l - cost 
a 



u 



(4.64) 



Now we consider the remaining integral 
1 a 



A 



d9 (1 + cos 6) exp 



2u2 



(l~cos6l) 



erfc 



V2(l- cos 61) 



(4.65) 



First we define (f) ~ 9/2, so that d9 = 2d(t). Thus l + cos6' ~ 2cos^ and similarly 1 — cos6' ~ 2sin^ 
In addition, since < 6* < tt, it follows that < <j) < tt/2, and hence both cos0 and sine/) are positive 
quantities. With this substitution we have 



A 



4 a r'^ 



cos" (p exp I — ^ sm (p eric — i 
/ V a 



(4.66) 



Using Eq. (2.31), we see that 



A = 



4 a r/\ 1 ( 2u 

A oo 1 / o \ " /■•"■/2 

4 a 1 ( 2u \ / ' 



sm ( 



cos sm 



(4.67) 



The integral over (f) is a standard one, and we find that 



A = 



1 a \ - 

V ri=0 



r(H±i; 



r( 






1 


+ - 


1 a 








1 


1 

+ - 


8 u 


U 




-v/tt q;^ 



2u 
a 



r(n±i) 



n— 3 



2m 
a 

2m 

a 



(4.68) 



Substituting this back into Eq. (4.64), we see that 
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At this point we can invert this Laplace transform term by term, using the fact that the inverse 
Laplace transform of a^'^ is f^^^ /T [v). We finally arrive at 



(5'(0,t)) = —=ut 



^ ^ormrmr(n + 2) 



-2utf 



(4.70) 



The even terms in this sum reduce to 



n=0 



r(M)r(M)r(2n + 2) 



{-2utY 



T{n + 2) 



Aut 

,^or(« + l)r(n + Dr(n + i)r(n + |) 

256 ^ 3 5 7 2 2 



(4.71) 



where in the third step we have used Eq. ( ^.30 ) and in the final step we have used Eq. (2.35) 
Similarly, the odd terms in this sum reduce to 



•S'n.odd — —j=ut 



- ^or(M)r(M)r(2n + 3) 



{-2ut) 



2n+l 



r(n + |) 



-4E 



^„r(n + 3)r(n + 4)r(n+|)r(n + 2 



^^r(n + 2)r(n + 3)r(n + i)r(n + i) 



(4.72) 



and so we have finally that 



(^(0,i)) = l 



256 
457r3/2 



ut 1F3 2; 



3 5 7 
2' 2' 2' 



ii^3 ( I i 2,3:^2^2 



(4.73) 



Figure 4.1 shows a plot of this result as a function of time. Looking back to Fig. 1.3, we see that 
this signal saturates too quickly to provide a good fit to the experimental data. In fact, looking ahead 



to Fig. 7.3, we see that it reaches saturation too quickly to favorably compare with the numerical 
simulation results for m = at resonance. 




Figure 4.1: A plot of the averaged signal as a function of time foi u = v and A = in the Hubbard approximation. 
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4.4 The coherent potential approximation 



Instead of making the Hubbard substitution of Eq. (4.33), we can replace 1/fi by its average in 
Eq. ( ^4.12b| ), obtaining 

f,^a + f2uii(jy (4.74) 



1=1 



This corresponds to a coherent potential approximation (CPA). 



4.4.1 Computation of (1//) 



We have from Eq. (4.74) that 



1 



1 



We again note that {1/ fk) is independent of k, so we drop the subscript and we see that 



(4.75) 



1 



The averaging procedure of Section 2.2 is by now second nature, and we see that 



(4.76) 



df3 exp 



13 (a+U^ 



dp exp I —/3a ~ ^y'^ \ / 



1 y/n u 



rfc 



(4.77) 



This is an implicit equation for (1//), which can be rearranged to give 



a 



//I 

2 x/^V\/ 



exp 



u 

4a 



erfc 



(4.78) 



Eq. (4.78) can then be solved numerically for (1//). 



4.4.2 Computation of (l//o) 



Just as we saw when considering the Hubbard approximation, the modification of Eq. ( 4.12b| ) to 
Eq. (4.74) makes the latter equation inconsistent with Eq. (4.12a). It is possible to use Eq. ( 4.78| ) 
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coupled with Eq. (4.12a) to obtain an expression for (l//o), but it does not reduce to Eq. (4.7S 



when we take the hmit A = and u = v. In this case, however, we can modify Eq. (4.12a) to 



/o = a + i A + V' 



so that it and Eq. (4.74) are again consistent 



Proceeding as we did in going from Eq. (4.76) to 4.77, we have 



and hence 



1 



dp exp 
d(3 exp 



/3 a + i A + 



-/3(a + ^A) - 



'TT 



V 



1 

a + iA 2 (a + iA)3/2 V\7 
■2/1 



exp 



V 



4 (a + i A) \ / 



erfc 



2Va + i A V \/ 



(4.79) 



(4.80) 



(4.81) 



4.5 Comparison of the Cayley, Hubbard, and CPA approxi- 
mations 

We have examined these approximations and found that the resuhs that agree best with the simu- 



lations are those obtained using the CPA approximation. In fact, as we will see in Section 4^, the 
CPA results agree strikingly well with the simulations. This is somewhat surprising, since in the 



Cayley equations ( 4.12a ) and ( 4.12b ) the on-site Green's functions for the U problem are allowed 
to vary from site to site, whereas in Eqs. ( 4.33| ) and ( 4.74 ) they are not. However, the Cayley tree 
approximation is exact in a one-dimensional system with nearest neighbor interactions and becomes 
increasingly inaccurate as the dimensionality of the problem and the range of the interaction in- 
crease. The CPA, on the other hand, is a mean field approximation, and it becomes accurate for 
higher dimensional problems and long-range interactions. 
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We see in Fig. 4.2 plots of Re (ao (Ajw)) for u = and several values of A. The Cayley and 
CPA approximations coincide in this case. The exact result for u — Q can actually be obtained 
analytically by the methods for computing ( S (A, a) \ used in Ref. [Rsl and described in detail in 



Chapter H. We reproduce here the result of Eq. (2.28), which is 



(ao (A, a)) 



a + iA 2 (a + iA)^/^V« 



exp 



4a (a + lA) 



erfc 



2y/a {a + iA) 



(4.82) 



As one would expect, this agrees with Eq. ( 4.32| ) for Q (a) — 1/y/a and with Eq. ( 4.81 ) for (1//) = 
1/a. 

We note the "gap" that is present in each of the curves, and we also note that the width of the 
gap increases with A. This effect can be seen analytically from Eq. ( 4.82| ). We start by converting 
from Laplace to Fourier space by making the substitution a = Lu/i, so that we have 



(ao(A,^)) = 



uj- A 



{lu - Af^ ^ 



exp 



4w (w - A) 



erfc 



2^uj {w- A) 



Using Eq. (2.31), we see that 



(4.83) 



Re (ao (A, w)) = Re 



2 (u.-Af'\/^ 



exp 



4w(w- A) 



(4.84) 



which reduces to 



Re (ao(A,w)) 



exp 



Au) (w - A) 



(4.85) 



if > A, and zero otherwise. This explains the gaps seen in the A — v and A = 2v curves. For 



A = 0, Eq. (4.85) becomes 



Re(ao(A,^)) = — -exp(- — 



(4.86) 



which produces the quasigap around a; = seen in the resonance curve. 



In Fig. 4.3, which presents plots of Re (ao (A, t)) for u = 0.25u and several values of A, we begin 
to see differences between the CPA and Cayley approximations. We see that the CPA results retain 
the gap structure of the u = case somewhat, while the only vestiges of the gap remaining in the 



Chapter 4: Sparse ss' — + pp' Processes with U Interaction 



74 



Cayley results are the dimple in the A = curve at lj/v — and a barely perceptible bend in 
the A — V curve at w/w = 0. The ofF-resonant CPA results are also noticeably asymmetric around 
w = A. 



In Fig. 4.4 we see plots for u = v, and many of the trends of Fig. 4.3 continue. The ofF-resonant 
CPA results continue to retain some memory of the gap and are still asymmetric around a; = A. 
The gap structure is no longer evident in the on-resonance curve. For A = the result for the 
Hubbard approximation is also shown. 



We see that the CPA and Cayley results look the most similar in Fig. 4.5, where u = 4w. However, 
upon careful examination one notes that the CPA results are still slightly asymmetric about a; = A. 
The peaks are also slightly wider for the CPA results. 

4.6 Plots 

4.6.1 The averaged amplitude (ao) as a function of time 



Figs. 4.6 and 4.7 display the numerical simulations of Re (ao (A, t)) and Im (ao (A, t)), respectively, 
for u/v — Q, 0.25, 1, and 4. We recall that (ao (A, t)) is the probability amplitude that the ss' — > pp' 
has not occurred during the time t. The striking feature of these graphs is that, off-resonance, 
oscillations with an approximate period A persist for a long time, not only when the exciton does 
not propagate {u — Q), but even for u = 4w. 

The graphs for u — Q can be obtained analytically by the methods for computing {S (A, t)) used 
in Ref. and described in detail in Chapter ||. In fact, we have obtained (ap (A,t)) by expanding 



Eq. (4.82) (or, equivalently, Eq. (2.2S)) in powers of 1/a and inverting the Laplace transform term 
by term. As in Chapter H, various formulas involving special functions can also be obtained. The 
exact analytical results in this case (not shown here) agree with the numerical simulations about as 



well as in Fig. 7.1. The initial dependence on time is given by 



(ao (A, i)) « 1 - ^vt - iAt (4.87) 

Looking at the graphs, it seems that the initial decrease of Re (ao (A, t)) depends on u, but this 
perception only indicates that the range of validity of the linear approximation Eq. ( |4.87 ) is small. 



On the other hand, the graphs of 1 — ( |ao|^ ) as a function of time in Figs. [7T1- |7.5| look linear and 
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Figure 4.2: Plots of Re (Reao {A,u>)) for u = and, from left to right, A = 0, and 2v. The A = v curve has been 
shifted vertically upward by 1, and the A = 2v curve by 2. The results are exact in this case. 
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Figure 4.3: Plots of Re (ao (A, a;)) for u = 0.25v and, from left to right, A = 0, v, and 2v. The A = v curve has been 
shifted vertically upward by 1, and the A = 2v curve by 2. 
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Figure 4.4: Plots of Re (ao (A,a))) for u = v and, from left to right, A = 0, v, and 2v. The A = curve has been 
shifted vertically upward by 1, and the A = 2v curve by 2. 




Figure 4.5: Plots of Re (ao (A, a;)) for u = 4v and, from left to right, A = 0, v, and 2v. The A = v curve has been 
shifted vertically upward by 1, and the A = 2v curve by 2. 
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independent of u over a wide range of vt. Note also that |(ao (A, w 1 — ^/^Tvt decreases initially 
twice as fast as (\aQ (A,i)|^^ w 1 — (y^/2) vi. Of course, it must be true that the mod square of 
the average is smaller than or equal to the average of the square, but the above result shows that 
|(ao (A, t))^ is not even a fair approximation to ( |ao (A, t)^ 



4.6.2 The averaged amplitude (ao) as a function of frequency 



Figs. 4.g- 4.11 present graphs of Re (ao (A, uj)) that correspond to the graphs of (oq (A, t)) in Figs. 4.6 
and |4.7| . For u — v and A = 0, Re(ao(A,a;)) is the exciton spectral density (times tt); more 
generally, it is an excitation spectral density for the resonant process ss' pp' . Each figure gives a 
comparison of the numerical simulation with the CPA approximation, which in particular for u = Q 
reduces correctly to the real part of the exact result of Eq. ( 4.82| ), with a — uo/i. We compare spectral 



densities, rather than time graphs, in part because obtaining (ao (A, t)) for the CPA approximation 
is laborious, and in part because the comparison of spectral densities is instructive, as we now 
discuss. 



In the notation of Section 7.3, the spectral density is obtained numerically from Eq. (7.7) as 



Re ao (A, lo + is) = Im V (4.88) 

in the limits e — > and N oo. In practice it is difficult to go beyond N — 10^, and one must 
choose £ to be fairly large to obtain a smooth plot for the spectral density where the eigenvalues are 
sparse, even after averaging over 10^ realizations. On the other hand, the distribution of eigenvalues 
for a random dipolar gas is highly structured near lu ~ 0, and if e is too large this structure will be 
missed. Rather than spending much time to go to large N and refine the numerical technique, we 
set N — 100, plot the results for e/v — 0.1 and 0.01, and let the eye extrapolate to the correct limit. 
It is remarkable that the CPA approximation agrees best with the smallest value of e/v, except 
for wiggles that are usually in regions of sparse eigenvalues. This is exactly how the exact spectral 
density is expected to compare with a calculation for finite N. A finite e can also be viewed as a 
simulation of the decoherence introduced by atomic motions, or of losses to other channels, such 
as the decay of the excited atomic states due to blackbody radiation as described in Chapter 5 of 
Ref. Thus, e is almost equivalent to the 7 of the Lorentzian broadening model of Ref. |l8|, the 



difference being that in Ref. ||l^ i"f is added to lu for the Fourier transform of Eq. ( 2.1b ) only, while 
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Figure 4.6: Plots of the real part of (oo (A, t)) for u = 0, u = 0.251;, u = v, and u = Av. 
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Figure 4.7: Plots of the imaginary part of (ao for u = 



= 0.25v, u = V, and u 



= 4v. 
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here ie is added to every w that appears. As e increases, the highly structured spectral densities of 
the ideal frozen gas become increasingly similar to Lorentzian lines. Because of unavoidable losses 
and incoherences, the experimental spectral densities will resemble the results of the simulations for 
some finite e. 



In Fig. 4.?, where results are presented for m = 0, the spectral density is quite different from a 



Lorentzian. At A = there is a quasi-gap around lj — Q. As A increases the gap extends from 



ix) = to = A, as was described in Section 4.5 



In Fig. 4.£ are presented results for u — 0.25w. We see that the sharp features of the u = 



spectrum are rounded off and the gap is beginning to be filled, but is still recognizable. 



Results for u = v are presented in Fig. 4.10. We note that the spectrum at A = has the 



appearance of a single, asymmetric line. At A = u and A = 2w a shoulder on the left of the line is 



the remnant of the gap of Fig. 4.8 



Fig. 4.11 contains results for u = Av. We see that there is a single line for any value of A, 
which becomes increasingly sharp as u increases. The lineshape is not Lorentzian, but even a small 
additional broadening will make it nearly so, as indicated by the calculation for ejv — 0.01. This 



sharp line leads to the damped sinusoidal oscillations seen in the u = 4i; graphs of Figs. 4.6 and 4.7. 
The general phenomenon of line narrowing can be discussed using the simple Lorentzian model for 
the exciton band, in which Eq. (lb) of Ref. |l^ is replaced by 

ick = Vfcao - i7Cfc. (4.89) 



Solving the equations in this limit, one finds a pole at 



' 7^ + A2 



(4.90) 



which gives a decreasing width as u increases, since 7 is proportional to u. One cannot simply 
replace 'Ylik average, which does not even exist. Instead, one can Fourier transform the 

average to find the time dependence 



exp (-iAt - \/leqt) , 



(4.91) 
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where 

= 7^2^ (4-92) 

and V has been defined in Eq. (2.12). A similar argument was applied in Ref. to find the time 

|2\ 



dependence of \ \an\ 

From the comparison with the simulations, it is apparent that a shortcoming of the CPA approx- 
imation is that it gives a symmetric spectral density at resonance (A = 0), while the correct spectral 
density is slightly asymmetric, except for u — 0. Correspondingly, (ao (A,t)) is not purely real at 



resonance, as can be seen in Fig. For A 7^ 0, the spectral density at w < A is underestimated 
by the CPA approximation. 

Clearly, the asymmetry of the spectral density at A = comes from the terms with I ^ m 



in the denominator of Eq. (4.4), which are neglected in the CPA approximation. It is surprising 
that these terms have such small effect on (ao (A,(jj)) at A = 0. At finite A, the underestimate of 
Re (ao (A, lo)) for w < A is probably due to the neglect of these same terms. We have already shown 
how to include these terms in the CPA formulation |jll|, but in practice, when the simple CPA is 
not adequate one can just as well use the numerical simulations. 



4.7 Summary 

In this chapter we have extended the simple model of Chapter |^ by allowing the sp ps process to 
take place. Although this complication no longer allows a direct analytical solution of the problem, 
we consider three approximations, each of which allows us to make some progress. We compare 
these approximations with the results of numerical simulations, which are carried out for a finite 
number of atoms by the methods described later in Chapter ^ of this thesis. 

The results of the Cayley tree approximation are not found to give the best results for the 
problem we are considering here, but they are still extremely worthwhile results. They extend the 
Cayley tree results of Rcf. ||l|], which are used extensively, for example, in the study of dipolar 



liquids H H, |TH|, 




The results of the Hubbard approximation are somewhat disappointing, since in this case we can 
obtain an explicit expression for the exciton spectrum (1//;) but we are not able to obtain a consistent 
result when the effects of the ss' —^ pp' process are included. However, in this approximation we 
can obtain the average of |ao|^ as a function of time when the V process is absent, and perhaps in 
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Figure 4.8: A plot of the real part of (oo {A,oj)) for u = and, from left to right, A = 0, v, and 2v. To aid the eye 
in separating the curves, the A = curve has been vertically offset by 2 and the A = v curve by 1. 
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Figure 4.9: A plot of the real part of (ao {A,u>)) for u = 0.25v and, from left to right, A = 0, v, and 
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Figure 4.10: A plot of the real part of (oo {A,u>)) for u = v and, from left to right, A = 0, v, and 
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Figure 4.11: A plot of the real part of (ao (A, ui)) for u = 4v and, from left to right, A = 0, v, and 
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the future this will lead to a more general technique that will allow this quantity to be computed in 
the other approximations as well. 

The CPA approximation is found to agree best with the results of the numerical simulations, 
and we are able to carry our computations in this approximation just as far as in the Cayley tree 
approximation. It seems at first strange that the CPA approximation gives better results than 



the Cayley, since in the Cayley equations (4.12a) and (4.12b) the on-site Green's functions for 



the U problem are allowed to vary from site to site, whereas in Eqs. ( 4.33 ) and ( 4.74 ) they are 
not. However, the Cayley tree approximation is exact in a one-dimensional system with nearest 
neighbor interactions and becomes increasingly inaccurate as the dimensionality of the problem and 
the range of the interaction increase. The CPA, on the other hand, is a mean field approximation, 
and it becomes accurate for higher dimensional problems and long-range interactions. It may also 



happen that, when the off-diagonal terms of Eqs. (4.4) and (4.11) are included, the Cayley is found 
to be the better approximation. 



Chapter 5 

Analogy With Spin Glasses 



5.1 Introduction 



In this chapter, we examine a correspondence between the atomic system we are seeking to model 



and spin glasses. We start in Section by considering the rubidium experiments of Anderson et 
al. H and Lowell et al. |Q. In Section 5.2.1 we derive in a rather straightforward manner a 



Hamiltonian that describes these experiments as a system of two interpenetrating spin glasses. We 



also determine the corresponding equations of motion. In Section 5.2.2 we illustrate how the results 



of Section 5.2.1 reduce to those of Chapter M in the appropriate limit. We derive a new, and seemingly 



more cumbersome Hamiltonian in Section 5.2.3 that also describes the rubidium experiments and 



is entirely equivalent to the Hamiltonian of Section 5.2.1. The method used to generate the new 



Hamiltonian is more generally applicable, however, as we demonstrate in Section 5.3 by using the 
new method to construct a spin Hamiltonian that describes the cesium experiments of Mourachko 
et al. Bs, H. Finally, we conclude in Section 5.4. 



5.2 The rubidium system as a spin glass 
5.2.1 Spin glass Hamiltonian and equations of motion 

We consider N atoms at positions which are initially in state s, and N' atoms at positions r^/ 
which are initially in state s' . We represent the s and p states at with the down and up states of 
an effective spin dfe , respectively, and similarly represent the s' and p' states at r^/ with a spin . 
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The Hamiltonian that describes this system with the ss' pp' and sp — s- ps processes mediated by 
dipole-dipole interaction potentials V and U, respectively, is then 

N N' 
k=l k' = l 

E ^fcfc' [^fc 4 + '^fc'] + E Uki<jt<jr^ (5-1) 



where 



and 



Vkk' = - — — — (5.2) 
rfc - rfc'l 



Uki = , ,3 , (5.3) 

ii I ^ k and Uki = otherwise]^ This is precisely the Hamiltonian for two interpenetrating spin 
glasses that interact via dipole-dipole potentials. 

For the experiments under consideration /i « 4/i', and the sp ps coupling p^, which leads to 
"spin diffusion," is larger than the ss' —>■ pp' coupling /i/i' that is responsible for the resonant energy 
transfer. 

We can write down basic equations of motion such as 



fc=i 



k' li^k 



and 

d ^ 

'Jt''k'-2T.^k'k{<yicT+-<^k<^k'), (5-6) 

k=l 

where A is the detuning from resonance.^ From these we can obtain 
1 

kl' kl 

- J2 ^k'kUkicyl (44 + ^fc^'^r) - A E ^k-k K'4 + ^k'^k) ■ (5-7) 

kA^k k 
^Of course, the potentials of Eqs. ( ^.2| ) and (5.3) are just the simplest examples of dipole-dipole interaction 



potentials. We are free to choose more complicated potentials if we so desire. 

^ To o btain these expressions we take the equations of motion that come about naturally from the Hamiltonian in 
Eq. ( [3.l| ) and introduce a phase exp [— i (cp — es)t\ into the definitions of cr^ and c^, . 
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These differential equations are rattier difficult to work; with, but we have already discussed one 
possible approximation scheme in Ref . p8] . 



5.2.2 Reduction to the sparse limit 



We now show explicitly how the general spin- variable formalism of Section 5.2.1 relates to the 



coefficients uq and of Eqs. ( 4.1a ) and ( 4.1b ) when there is only a single atom that can be in the 
s', p' pair of states (i.e., when only a single spin variable of the type ak' is present.) The key to the 
correspondence is that in this case the initial state \i), which consists of all spins down, evolves to 



\t) = Uo(A,0 + Ec9(^'^)'^^< I*) 



(5.8) 



at time i, where the spin variables now denote time-independent Schrodinger operators. With the 
understanding that r^' =0, we have: 



(5.9) 



and 



{t\al\t)^2\ck (A,i)|^-1, 



It is then clear that the expectation value of Eq. (5.6) is equivalent to 



|ao (A, t)f = J2 [4 (A, t) ck (A, t) - 4 (A, t) ao (A, t) 



(5.10) 



(5.11) 



(5.12) 



where T4 is short for 14' fe. This equation is itself a consequence of Eq. ( 4.1a ). Similarly, the equations 
for {t |(7|| t) and |o'fc'0'fe^| t) are equivalent to the equations for |cfc (A,i)|^ and Ck (A,t)^ ao (A,t) 



that follow from Eqs. ( 4.1a ) and ( 4.1b ) 



5.2.3 Another spin glass Hamiltonian for the rubidium system 



In the rubidium experiments of Anderson et al. |6j and Lowell et al. there are the four states 
s, p, s', and p' . Therefore we can define a column vector of length four for each site, so that if an 
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atom at site A is in the s, p, s', or p' state, then the column vector corresponding to that site is 
(0,0,0,1)"^, (0,0,1,0)^, (0,1,0,0)"^, or (1,0,0,0)"^, respectively. The Hamiltonian for this system 
can then be written 



A Ai,A2 L " 

J2 f^AiA2 [«)a, K"p)a. 



+ 



Ai,A2 



where e, afp, and afi^, are the 4x4 matrices 










































V 










J 



(5.13) 



(5.14) 



and 













V 



a 










(5.15) 



(5.16) 



and the subscripts on these matrices indicate on which site they are to be applied. The quantities 
VaiA2 and Uk^k^ are defined to be zero if Ai = A2, and are the interaction potentials 



|rAi - rAj 



|3 ' 



and 



t^AiAa = 



M 

|rAi - rA2 



(5.17) 



(5.18) 
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otherwise. In the interest of completeness, we note that we could also include the s'p' — > p's' 
interaction by adding a term 



Ai,A2 



(5.19) 



to the Hamiltonian of Eq. (5.13) 



The Hamiltonian of Eq. (5.13) is entirely equivalent to the Hamiltonian of Eq. (5.1). The fact 



that Eq. (5.1) contains one set of spins corresponding to the primed sites and another set of spins 



corresponding to the unprimed sites is seen to arise from the fact that Eq. ( 5.13 ) can be viewed 
as containing two distinct types of operators. One type operates only on the s' and p' states and 
hence only on primed sites, while the other acts only on the s and p states and hence only on the 



unprimed sites. Eq. (5.13) appears more complicated than Eq. (5.1), but it actually embodies the 
general technique for constructing spin Hamiltonians for these systems, as we will see in Section f3.3|. 



5.3 Spin glass model for the cesium system 

In the cesium experiments of Mourachko et al. there are three states that contribute to 

the experimental results. These are the 23s, 23p, and the 24s states, which we refer to as the s, 
p, and s' states, respectively. The processes pp — > ss', ps — > sp, and ps' s'p can occur. The 
pp — > ss' process can be tuned in and out of resonance via the Stark shift by application of an 
electric field, while the ps — > sp and ps' s'p processes are always resonant. In the case of the 
cesium experiments, the ps — > sp and ps' —^ s'p processes are of comparable strength and both must 
be considered. 



Just as we did in Section 5.2, we define a column vector of length three for each site. If an atom 



at site A is in the s, p, or s' state, then the column vector corresponding to that site is (0, 0, 1) , 
(0, 1,0)"^, or (1,0,0)"^, respectively. The Hamiltonian for this system is then 



Cs 



E 



eA 



Ai,A2 



Ai,A2 



J2 ^AiA2 [K.)a^ ('^^Ja, +h.C. 
Ai.Aa 



(5.20) 
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where e, a^^, , and a^^, are the 3x3 matrices 



ep 




(5.21) 



(7 , = 

ps 



and 



V 



















(5.22) 



(5.23) 



and the interaction potentials VaiA2i Ua-^a^, and WaiA2 ^-re zero if Ai = A2 and otherwise are the 
potentials corresponding to the pp ss', ps' — > s'p, and ps — > sp processes, respectively. 

The equations of motion can be determined in this case as well, simply by computing commuta- 
tors. 



5.4 Summary 

We have seen in this chapter how the frozen dipolar gases studied by Anderson et al. ^ and 
Lowell et al. |^^, as well as Mourachko et al. |4|, are completely equivalent to systems of 
interpenetrating spin glasses. Casting the problem in this form has the advantage that it becomes 
possible to at least write down equations that go beyond the sparse model of Chapters and 0. 



Chapter 6 

The Effect of Spin 



6.1 Introduction 

Throughout most of Chapter 0, we used an interaction potential that varied as 1 jr^ , with no angular 



dependence. In Section 2^ we discussed how to deal with interaction potentials that have a more 
complicated dependence on r, depend on the angular variables, or both. In this chapter we determine 
exactly what the spin-dependent interatomic potentials are for the ss' ^ pp' and sp — s- ps processes 
in the experiment we are modeling ^, The interaction potential depends directly on the rrij 
values for the states involved, where nij is the projection of the total angular momentum in the 
direction of the applied electric field. For states with I = 0, nij is just the projection of the outer 
electron's spin. For states with 1 = 1, the spin-orbit coupling separates the states with j = 3/2 from 
the states with j = 1/2. Further, we assume that for j = 3/2, the applied electric field splits the 



states with \mj\ =3/2 from those with \mj\ = 1/2. As is illustrated in Fig. 1.1, this is the case for 
the 34p states of Rb under the conditions of the experiments in which we are interested. Then, for 
each of the states we consider, \mj\ is fixed and rrij can only take on the values ± \mj\, so that it is 
effectively equivalent to a spin of 1/2. When we refer to spin effects in this chapter we mean, more 
precisely, the effects due to rrij which are derived from the spin-orbit coupling. 



We start in Section x2 by examining the matrix elements of the dipole-dipole interaction in 
more detail, and taking into account the effect of the various atomic spin states. The techniques 
we use for doing this are inspired by Ref. p^ , which we found to be extremely helpful in that it 
pointed us in the right direction when we began our study of spin effects in the frozen system. In 
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Section 3.2.1 we recall the form of the dipole-dipole interaction potential and in Section 6.2.2 



begin the calculation of the dipolar matrix elements with the effects of the atomic states taken into 
consideration. This computation requires a change of basis of the atomic spin states, as explained in 
Section 6.2.3. In Section [6.2.4 the laborious procedure necessary to compute the individual matrix 
elements is explained, and the results are presented in tabular form. 



Next, in Section 6.3 we attempt to determine the effective interaction potential for the resonant 



and nonresonant processes. We start in Section 6.3.1 by determining the states and probability 
amplitudes of the system when the complications due to the presence of atomic spin are taken into 



account. We proceed in Section 6.3.2 to discuss the corresponding equations of motion. Wc find, 
as expected, that in general the quantity rrij is not conserved. Therefore the number of interacting 
states for a system of + 1 atoms is 2^+^, whereas it is only iV + 1 when the effects of spin 
are neglected. However, when states of \mj\ = 3/2 are involved it is still possible to obtain exact 
formulas as we did in Chapters, provided that the quantity V is replaced by an effective interaction. 



We refer to this case as the separable case, for reasons that will become apparent in Section 5.3.4. 



Section 6.3.3 concerns the properties of Dirac matrices, which are instrumental in our discussion 



of separability and its applicability to averaging in Section 6.3.4. In Section 6.3.5 we determine 
the effective interaction potential in the separable case. Finally, we summarize and conclude in 



Section 6.4. 



6.2 Matrix elements for the ss' pp' and sp ps processes 
6.2.1 The dipole-dipole interaction 

We consider two atoms, the first with its nucleus located at x„uc and outer electron located at Xdoc, 
and the second with nucleus located at x^^^, and outer electron located at x^j^^. Next we define 

^nuc -'^nuc (^■-^) 



^ — ^clcc ^nuc; (^■^) 

and 



dec 



(6.3) 
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The potential for the electrostatic interaction between the two atoms is then, in atomic units, 



C/(r,r',R) = - 



1 



R R + r' 



R + r' R-r 



(6.4) 



The first two terms on the right hand side represent the interaction between the two nuclei and the 
interaction between the two electrons, respectively. The last two terms represent the interactions 
between a nucleus and the electron that belongs to the other atom. 



As is shown on page 283 of Ref. pa and in more detail on page 260 of Ref. |55 , expanding 



Eq. (6.4) and keeping terms only to leading order in r/R and r' / R gives 



C/fr,r',R) 



r • r' - 3 R • r R • r' 



i?3 



(6.5) 



This is the well known result for the dipole-dipole potential between two atoms. 



It is important to note that, although Eq. (6.5) is the same potential used in the derivation of 
the van der Waals interaction, the effects we study are transient effects that take place after the 
atoms have been allowed to interact only a very short time. Therefore they are quite different from 
the steady-state van der Waals effect. 



6.2.2 Dipolar matrix elements 

We will later be averaging over the positions of the atoms, so we only want to determine the angular 
dependence of the matrix elements. Thus the 1/R^ factor can be ignored until averaging, and we 
need only consider matrix elements of the type 



r • r' - 3 (r • r) (r' • R.) s'^ 



(6.6) 



• r' - 3 (r • r) (r' • R.) . (6.7) 

Here we note that when we write the states \a) \b) and (c| {d\ we really mean \a) (E) \b) and (c| (g) (c?|. 



respectively. In Section 6.2.4 , where we deal with products of a great many states and the tensor 
products are more urgently called for, we will write the tensor product symbols out explicitly. For 
now, however, since we are only dealing with products of two states at a time, we omit them to save 
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space and make the notation a little simpler. 

We now recall that, in the experiments we are modeling, s, s', p, and p' refer to the 25s, 33s, 
24pi/2, and 34p3/2 states, respectively. We also recall that there are really two distinct resonances 
in the ss' — > pp' interaction, corresponding to whether \mj\ — 1/2 or \mj\ —3/2 for the p' atom. 
The complication for the interaction potential enters through the p and p' states. Clearly the matrix 



element in Eq. (6.7) will be different for the case where the p state on the left has mj — 1/2 and 
the one on the right has nij — —1/2 than for the case where both p states have rrij — 1/2, because 
the p state with rrij — 1/2 will involve a different mixture of the p^, Py, and p^ states than the state 
where rrij — —1/2. The matrix element in Eq. ( ^.(]| ) will similarly depend on the rrij values of the 
p and p' states, although there is some simplification because we are able to consider the p' states 
with \mj\ — 1/2 and those with \mj \ = 3/2 separately. 
For simplicity, we define 

p = r • r' - 3 (r • r) (r' • r) , (6.8) 

and we note that 

p — xx + yy + zz — 3 + yRy + zR^ {x R^ + y Ry + z'R^ , (6-9) 

where Rx — sin cos 0, Ry — sin sine/), and R^ — cos 9. 
Consider for a moment the matrix element 

{Pa\{p't \XX'\ S')| S) = {pa \x\ s) {p[ \x'\ s') , (6.10) 

where a and b can be x, y, or z. We can apply the usual parity argument to determine that {pa \x\ s) 
is nonzero only if a = x, and in that case is equal to the dipole moment fi. Similarly, (pj, \x'\ s') is 
nonzero only ii b = x, and in that case is equal to the dipole moment n' . Thus we have 

{pa\{p',\xx'\s')\s) = {pa\x\s){p',\x'\s') 

= PfJ-'SaxSbx- (6-11) 
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Analogously, 

{pb\{s\xx'\pa)\s) = {pb\x\s) {s\x'\pa) 

= IJ'^SaxSbx- (6.12) 

In fact, the matrix elements we consider are a little more complicated than this. The atoms have 
a total spin J = L + S. The potential p changes L, but does not at all affect S and therefore cannot 
change s or . Therefore the matrix elements we really consider are of the form 

{Pa,a\{Pb,f3\xx'\s',j)\s,£) = {pa,a\x\s,e) {p'i„(3\x'\s',j) 

= PP''SaxSbxSaeSi3j, (6.13) 

where, for example, \s,a) is short for \s,mg = a). Analogously, 

{pa,a\{s,P\xx'\pb,j)\s,e) ^ {pa,a\x\s,e) {s,(3\x'\pb,'^) 

= p'^SaxSbxSaeSpj. (6.14) 



Using Eqs. ( |6.13 ) and ( |6.14 ), it is trivial to see that 



{Pa,a\{p^,(3\p\s\j)\s,e) ^ pp! {Sab-'iRaRbj SaeSfi-y, (6.15) 

and 

{Pa,a\ (S, (3 \p\ Pb, 7) I S, e) = fl^ (^Sab ~ SRaRb^ SaeSi3-r- (6.16) 

6.2.3 A change of basis 

We can easily write out the states s, s' , p, and p' in terms of states of the form |j, Z, urij). For exam- 
ple, the state s denotes the pair of states |j = 1/2, 1 = 0, rrij — 1/2) and \j ~ 1/2, 1 = 0, mj = —1/2), 
while the state s' denotes the pair of states \j' = 1/2,1' = 0, rrij = 1/2) and \j' = 1/2,1' = 0, nij = —1/2). 
Similarly, p denotes the pair of states \j — 1/2, 1 — 1, rrij = 1/2) and \j ~ 1/2, 1 ~ 1, mj = —1/2). Re- 
calling that the p' states with \mj\ — 3/2 are treated separately from those with \mj\ = 1/2, we see 
that p' denotes either the pair of states \ j' = 3/2, 1' — 1, mj ~ 3/2) and \j' — 3/2, 1' — 1, rrij = —3/2), 
or the pair of states \j' ~ 3/2,1' ~ l,mj ~ 1/2) and \j' = 3/2,1' = l,mj = —1/2). In order to use 
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State 


Decomposition 




\j = h^^j = 1) 
i = ^,mj = — 1) 


-72 






[2 

V 3 


P^'5)-75(|P-'~5)+*l 

h'-^) + 7i(k'5)-« 


\Pv^l)) 




\j = l^rrij = i) 

J = = -1) 


i 

V3 


P-'-5)-73(l^'-'5)-^ 


Py^)) 



Table 6.1; The decomposition of all p and p' states into linear combinations of states of the form \pa,ms). 



Eqs. ( 6.15 ) and ( |6.16|) , we have to write the p and p' states as linear combinations of states of the 
form \pa,ms). This is easily done by consulting a table of Clebsch-Gordan coefficients, as can be 
found on pages 76-77 of Ref. |Q. Using such a table, we determine the decompositions in Table pT] . 



6.2.4 The ss' pp' and sp ps matrix elements 

With the decomposition of Table |6.l| and the rules of Eqs. ( 3.15 ) and ( 6.16| ) in hand, it is now 



possible to compute the matrix elements we need. I will first show how to compute one matrix 
element, as an example, then present all the matrix elements in table form. 
Consider the matrix element 



p,mj 



P 

mm' 



(6.17) 



We have from Table 6.1 that 



P ,rn^ 



V2 ( 



Px^- 



(6.18) 



and 



p,mj 



1 
73 



V3 ( 



Px,^ ) - I 



Dropping the occurrences of "rrij that appear to save space, we see that 



(6.19) 



P, 



1 






1 




2 






iV3 


Pz, 



1 

L 71 



1 

71 



Px,^)+i 



Px, 

1 

' 2 



Py^ 
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p/ (mm') 








' 2) 2) 




-^1 cose sin 


-^f sin^ 6'e2*'^ 








{p,-h 


(pM 




-^(l_3cos20) 


y|cos6'sin6'e*'^ 








{p^\ {p'^-'i 










cos sin 


^(l-3cos2^?) 


{p^~\ 














-y|cos6'sin6le-*'^ 



Table 6.2: The matrix elements {p,mj \(p',mj \p/ {lJ.fi')\ s',mj)\ s,mj) for the ss' —* pp' process when \mj\ = 3/2 for 
the p' atom. 













+ ^(py^^ 


)] 


X 



^-2 



^-2 



+ ^{Py^2 



HPx-^2 



(6.20) 



Now we simply apply Eq. (6.15) to each matrix element corresponding to a term in the sum of 



Eq. (6.20). For example, the first term in the sum corresponds to the matrix element 



1 




P 




2 


{p.,l 




4) 



(6.21) 



which according to Eq. ( 3.15| ) is equal to 



- sm W cos y cos ( 



(6.22) 



Computing the other terms similarly and summing them, we determine one of the matrix elements 
for the ss' — > pp' process. 

The process is straightforward, albeit laborious, and one obtains the matrix elements presented 



in the Tables 6.2 - 6.4 



6.3 Determining the potential 
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p/ W) 




2) ^' 2) 


2) 2) 


2) ^' 2) 


{p^ h {p'-^ \ 


-^(l-3cos2 6) 


V2cos6'sin6'e^'^ 


-^cos6'sin6'e*'^ 


— i= sin^ Oe^^'t' 


{p^-h 


{p'l 




V2cos6'sin0e-''^ 


^(l-3cos2 e) 




^ cos 6' sin ^e*"^ 


{p^-2 {p'^~\ 




^cos6'sin6'e-''^ 




_^^(l_3cos2 0) 


72 cos 6* sin 


{p^-\ 


{p'^-\ 




sin^ 616-2*"^ 


-^cos6isin6ie-^'^ 




^(l-3cos2 e) 



Table 6.3: The matrix elements {p,mj \{p',mj \p/ s',mj)\ s,mj) for the ss' — » pp' process when \mj\ = 1/2 for 

the p' atom. 













{pI 




i(l-3cos2 0) 


- cos e sin ^e-*"* 


- cos 6* sin 


- sin^ 6'e-2»<^ 


{pA ('^'-I 




— cos e sin ^^e*"^ 


(1 - 3cos^ e) 


(1 -3cos^6l) 


cos 61 sin ee-"f' 


{p^-i 






— cos e sin 06**^ 


(1 - 3cos^ 0) 


-i (l-3cos^i9) 


cos sin 


{p^-l 






- sin^ ^e^'"^ 


cos 9 sin 6e'"^ 


cos ^ sin ^e"^ 


^(l-3cos^ 61) 



Table 6.4: The matrix elements (p, rrij |(s, |p/ (/^^) | p, rrij) | s, raj) for the sp — > ps process. 



6.3.1 States and probability amplitudes with spin included 

In the initial state, the system is labeled by an index |s',mj) representing the state of the s' atom 
and a set of N indices {|s,mj)} representing the states of the N s atoms. In the absence of an 
applied magnetic field, the initial values of these indices are random. The state corresponding to 

the configuration where the primed atom is in the s' state with = ao and the unprimed atoms 
are all in the s state with {mj} = {ao; ai, . . . , ajv} is 

|s', ao) ®\s,ax)®---® |s, ajv) . (6.23) 

Of course, because the s and s' states have total spin 1/2, all are required to have the value 1/2 
or —1/2. The probability amplitude corresponding to this state is 

(i) = {(s, ajvl ® • • • ® (s, ai| (g) (s', ao|} \i3 (t)) , (6.24) 

where {t)) is the state of the system at time t. 

Because the system is allowed to make ss' — > pp' transitions, we must also consider states in 
which the primed atom is in the state p' and a single unprimed atom is in the state p. This state 
is labeled by the indices \p',mj) and \p,mj), together with the indices {\s,mj)} corresponding to 
the s atoms not involved in the ss' — > pp' transition. The state corresponding to the configuration 
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where the primed atom is in the p' state with rrij — Pq and all the unprimed atoms except for the 
fcth are in the s state while the fcth unprimed atom is in the p state with rrij — is 

\p', po) (g) |s, • • • (g) |s, (3k-i) \p, Pk) <E) \s, Pk+i) ®---®\s,Pn)- (6.25) 

The quantities /3i, . . . , Pk+i, ■ ■ ■ , Pn are restricted to the values ±1/2 because the s state has 
total spin 1/2. Because the p state also has total spin 1/2, Pk is similarly restricted to the values 
±1/2. The p' state has total spin 3/2, and the quantity Pq is restricted to the values ±3/2 or ±1/2, 
depending on whether we are considering the \mj\ — 3/2 or \mj\ — 1/2 resonance, respectively. 
Similarly to Eq. ( |6.24 ), the probability amplitude corresponding to the state in Eq. ( |6.25 ) is 



Ck:{^3o■,|3^,...M W = {{s, Pn\ ® ■ ■ ■ ^ (s, Pk+l\ » (p, Pk\ ® (s, Pk~l\ » ' ' ' ® {p , Po\} H W) , (6.26) 

where j"!/; (t)) is again the state of the system at time t. 

6.3.2 Equations of motion 

In the interest of brevity, we denote the subscript lists 

{ao; aw}, (6.27) 

and 

{7o;/3i,...,/3fc_i,7fc,/3fc+i,...,/5jv} (6.28) 

by A and i3fc(7o,7fc)j respectively. We label the 0th and fcth components of Sfc(7o,7fc) with a 
different symbol than the rest because we will sum over Bk (7o,7fc), and it will turn out that the 
Bk (7o,7fe) term in the sum corresponds to the s' atom at the origin and the s atom at position Vk 
making the ss' pp' transition. Thus the use of a different symbol serves as a reminder of which 
atoms are involved in the transition. 

The equations of motion for the a and c coefficients of Eqs. (6.24) and Eqs. ( |6.26 ) in Fourier 
space are then 

N N 

{u;-A)aA = *n'^o„.o<;'+E E ^^*^k(7o.7.)^^^^'=(^o,7.). (6-29a) 

m=0 k=lBk(lo,lk) 
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^Cfe;B,(70,7.) = Z]^'=^b7(70,7.),>1'"-4'' 
A' 



(fc) 



(6.29b) 



where aJn is the initial rrij value for the mth atom and 



(6.30) 



is just the radial part of the interaction potential. Here we have defined A' = {aQ;a[, . . . ,Q;^} 



to be another set of indices subject to the same constraints as A. The matrix M 



is a 



2^+1 X 2^+^ matrix that is constructed from the 4x4 matrix formed by Table x2 or |6.3| , depending 
on whether we are considering the case where the p' atom has \mj\ = 3/2 or \mj \ = 1/2. Explicitly, 
if the set of 

^ ' ^ ' ' ^ ' ^ (6.31) 



i{/3,/3'},{a,a'} {Ok,<t>k) — (p^ P 





p 











are the elements belonging to the 4x4 matrix that corresponds to the value of \mj\ that we want 
for the p' state, then 



M 



N 

l=l,l^k 



(6.32) 



We see immediately from Eqs. (6.32) and (3.31) that 



M 



(k) 

-4,-Bfc(7o,7fc) 



M 



(fe) 

-Bfc(7o,7fc),^ 



(6.33) 



Multiplying both sides of Eq. ( |6.29b| ) by VkM^^^^^^^^^^^^ and summing over k and Bk (70, 7fc), we 
see that 

N N 

E ^fe<k(70,7.)^feA(70.7.)-EE^M'i 
fc=lBfc(70,7fc) A' 

where 

-Bfc (7o,7fc) 



A, A' A,Bk{-yo,lk) -Bfc(7o,7fc).^ 



Using Eq. ( 6.33 ) we see that, in terms of matrices, Eq. ( 6.35 ) is equivalent to 



(6.34) 



(6.35) 



(6.36) 
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Plugging the result of Eq. (6.34) into Eq. (6.29a), we see that 



uj (uj — A) aA — ito 



m=0 fc=l A' 



(6.37) 



For iV + 1 atoms, this represents 2^+^ equations. 

The quantity ua actually depends on the two sets of indices A and yl^*' = |cto'^ • • ■ ; |j 
as is clear from the first term on the right hand side of Eq. ( |6.37 ) . Therefore we can consider a to 
be a matrix, and we can write 



N 



UJ {uj — A) a = iuj + 



(6.38) 



fc=i 



At this point we solve this equation formally and switch to Laplace space by making the substitution 
to — ia. We find 



aia + ^A) + Ek=lViP<^''^ 

a dfi exp [—(3a (a + iA)] exp 
Jo 



N 



k=l 



(6.39) 



Because Eq. (6.29b) is coupled to Eq. (6.29a), Cf^.Bkiio-ik) ^^^^ depends on the set of indices A'^ 



i) 



Thus Cfc can also be considered a matrix, and we can rewrite Eqs. (6.29a) and (6.29b) in the matrix 
form 



{to — A) a — i + 



N 



t 



fc=i 



U)Ck 



(6.40a) 
(6.40b) 



We note that these equations correspond to Eqs. (2.1a) and ( 2.1b| ), and so we hope to obtain an 
expression for S{A,a) analogous to Eq. (pT7|). We have to exercise a little caution, however, since 
Eqs. (2.1a) and ( 2.1b[ ) are equations for scalar quantities while Eqs. ( ^.40a ) and (6.40b) are matrix 
equations. First consider Eq. ( |6.39| ) in the rational form we had before we introduced the integral 
over p. Because the matrix J2k=i ^k-P'^''^ Hermitian we can diagonalize it. If we let |A) be the 
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eigenstate of ^^.^^^ V'^2p(fc) ^[^]^ eigenvalue A, then 



(AHA) = 



a{a + iA) + ' 



(6.41) 



The quantity (A \a\ A) is just a scalar, so we can certainly invert this Laplace transform, take the 



mod square, and Laplace transform again just as we did in Section 2.2. In this way we find that 



5* (A, a) 



A = 



2a 



1 - 



a^ + A^ 



>;2 + A2 + 4A2 



(6.42) 



We could do this for any of the eigenvectors |A), and so it follows that 



5(A,a) . i- 



1 - 



a 



1 a 
2^ 



a2 + A2+4^;L,lfp(fc) 

2 , ^2 



2a 



■ / d/3 exp [-P{a^ + A^)] exp 

^0 



N 



(fc) 



fe=i 



(6.43) 



which looks very much like Eq. (2.7). 

To obtain results that can be compared with the experiments, we also want to average over initial 
spin configurations and sum over final spin configurations. Thus we want to compute the quantity 



1 



-TrS'(A,a) 



1 a' + A ' 1-°° 
2^ 



2a 



/ d(3 exp[-p{a^ + A^)] 
Jo 



1 



^Tr exp 



N 



-4/3^T42p(fe) 



k=l 



(6.44) 



6.3.3 Dirac matrices 



Before we go any further with Eq. (6.44), we first examine some properties of 4 x 4 matrices. It will 



turn out that this will aid us greatly in carrying out the average of Eq. (6.44). 

We know from pages 211-213 of Ref. Q that any 4x4 matrix can be written as a unique linear 
combination of the sixteen Dirac matrices. These matrices are of the form CTq, (8" cr/3, where a and /3 
range from zero to three. Here we define ctq, to be the 2x2 unit matrix if a = and the usual Pauli 
matrix otherwise. Thus, if we have any 4x4 matrix Q then we can write 



Q = ^ ^ <7a/3 o-Q cr/3, 
Q=0 (3=0 



(6.45) 
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Pss' — >pp' .3/2 




2) 2) 


' 2) 2) 










f sin" 6 


-i sin (261) e-"^ 










- -^'^ ^ 




-isin (26*) e*"^ 


[5 + 3 cos (261)] 








{^^. 












[5 + 3 cos (261)] 


^ sin (26*) e-*"^ 


(^^4 












i sin (26*) e"'' 


1 sin" 6 



Table 6.5: The P matrix for the ss' pp' process when \mj \ = 3/2. 



-Pss' — >pp' ,1/2 


5,2) s, 2) 


^'-^> ^.-^) 


s , 2) 2) 


S , 2) 2) 








^ [29 + 3 COS (261)] 


sin (261) e-*-^ 


-:|sin(26')e-** 


-|sin" 6*6- 




! \* ' 2 




-i sin (261) e"f' 


4 [5 + 3 cos (261)] 


^ [5 + 3 cos (261)] 


i sin (26*) e-*"^ 








-| sin (261) e'f 


i [5 + 3 cos (261)] 


^ [5 + 3 cos (261)] 


^ sin (26*) e-*"?' 








-fsin" 6le"*"^ 


i sin (261) e"f' 


i sin (261) e"^ 


^ [29 + 3 cos (261)] 



Table 6.6: The P matrix for the ss' — » pp' process when \mj \ = 1/2. 



for some coefficients qap- These coefficients are in fact easy to determine, since 



1 



Qap = ^Tr [Q (aa •Sicrp)] 



(6.46) 



Eq. (5.46) follows easily from the well-known properties of the Pauli matrices and the facts that 

{A(g)B)iC(g)D) = {AC)(g){BD), (6.47) 



and 



Tr {A(g)B) = (TtA) (TrB) 



(6.48) 



Let Mgg,^pp, 3, Afgs'^pp'.i; s-iid Msp^ps be the 4x4 matrices determined by Tables |6^ 
Then define 



M 



SS — >pp 



(6.49) 



and similarly for Pgs'^pp' i ^nd Psp^ps- We will see in Section 6.3.4 that these three matrices are 
very important when it comes to averaging Eq. ( |6.44| ). Therefore we present the three P matrices 



in Tables 6.5-5.7 



We will also see in Section 6.3.4 that in order to obtain the average of Eq. ( 6.44 ) we will need 
the decompositions of the P matrices into linear combinations of the Dirac matrices. We have seen 
that we can use Eqs. ( |6.45 ) and ( 6.46| ) to determine these linear combinations, and in this way we 
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p 


P^9.) S,^} 






P.-k) s,-^) 




, 




i [7 - 3 cos (261)] 


-i sin (261) e-"^ 


-i sin (261) e-'* 


-|sin^ 9e-^"t' 








-isin (26*) e"f' 


4 [5 + 3 cos (261)] 


4 [5 + 3 cos (261)] 


i sin (261) e-^'f' 




{P^-k 




-|sin (26*) e^'f' 


4 [5 + 3 cos (261)] 


^ [5 + 3 cos (261)] 


1 sin (261) e-^'f' 








-|sin^ 6*6^'* 


i sin (261) e'* 


i sin (26*) e'f 


i [7 - 3 cos (261)] 



Table 6.7: The P matrix for the sp — > ps process. 



Dirac Matrix 


Projection 


CTo <i 


5 (7q 


^ [7 - 3 cos (261)] 


CTO 


5 CTi 





CTO 


5 (72 





(To ^ 


$ (73 





CTl 


5 (70 





CTl 


5 (71 





(Tl 


5 (72 





(Tl $ 


5 (73 





(72 $ 


5 (70 





(72 (< 


?1 (7l 





(72 1 


5 (72 





(72 (< 


5 (73 





73 <i 


5 (70 





73 <i 


5 (71 


— ^ sin {29) cos 


73 <i 


5 (72 


— ■^ sin (26*) smc/) 


73 <i 


5 (73 


i [1 - 3 cos (261)] 



Table 6.8: Projections of the 4x4 matrix P onto the Dirac matrices for the ss' — > pp' process when \mj \ = 3/2 for 
the p' atom. 



obtain Tables |6^^-|6lC . 



6.3.4 Separability 

At this point we would like to write 



exp 



= ]Jcxp 

fc 



-/3yfe2p(fc) 



(6.50) 



so that we can perform our usual averaging procedure as described in Section 2.3. Unfortunately, 



Eq. (6.5C) does not hold unless the commutator of P^'") with P^'^ is a c-number. We can use 



Tables 6.8-6.1C to check if this is the case for our three potentials. 

Because it is the simplest, we first consider the potential for the ss' — s- pp' process when \mj\ 
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Dirac Matrix 


Projection 


(To (To 


i [3 + cos (20)] 


(To (S) CTl 





(To (T2 





(To (T3 





(Ti (g) (To 





(Tl (El CTl 


1 

18 


5 + 3 cos (20) - 6 sin^ cos {2(j)) 




CTl (g) 0-2 


-isin^6'sin(20) 


(Tl (T3 


— ^ sin (20) cos^ 


(T2 (To 





(T2 (Tl 


-isin^0sin(20) 


(T2 (T2 


1 

IS 


5 + 3 cos (20) + 6 sin^ cos {2(j)) 




(T2 (T3 


— ^ sin (20) sin (f) 


(T3 (To 





(T3 (Tl 


— g sin (20) COS0 


(T3 (T2 


— ^ sin (20) sm(f> 


(T3 (T3 


[1-3 cos (20)] 



Table 6.9: Projections of the 4x4 matrix P onto the Dirac matrices for the ss' — » pp' process when \mj \ = 1/2 for 
the p' atom. 



Dune Maiiix 


Pivjcclion 


(To (To 


'2 


(To (Tl 





(TO (T2 





(To (T3 





(Tl (To 





(Tl (Tl 


1 

IS 


5 + 3 cos (20) - 6 COS {2(j>) sin"^ 




(Tl (T2 


-isin^0sin(2(^) 


(Tl (T3 


— :|sin(20) cos^ 


(T2 (To 





(T2 (Tl 


-gsin^0sin(2(^) 


(T2 (T2 


1 

IS 


5 + 3 cos (20) + 6 cos (2(^) sin'' 




(T2 (T3 


-isin (20) sinci) 


a:s ff(i 





(T3 (Tl 


— ■^ sin (20) coS(^ 


(T3 (T2 


— ^ sin (20) sin (j) 


(T3 (T3 


i [1-3 COS (20)] 



Table 6.10: Projections of the 4x4 matrix P onto the Dirac matrices for the sp ps process. 
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3/2. Using Table 3.S we see that 



p{k) ^ 



(fc) 
Pqq 



for some coefficients Pqq"* , p'^^^ , p'^2 : ^^'^ P33 that depend on 9k and 0^ . We also define 



.(1) 



(fc-i) (fc) (fe+i) 



(6.51) 



(6.52) 



Similarly, we have 



P 



(0 
P32 



(6.53) 



Because ctq is just the 2x2 identity matrix, it follows that {(7^°^ ® 4*-') } is just the 2^+1 X 2^+1 
identity matrix. Thus ^crl^"' <8) o'o'^''| will commute with anything. Now consider the product 



(6.54) 



Although the result does not depend on it, we assume for convenience that k < I. We also expand 



each of the terms in Eq. (6.54) and apply Eq. (6.47). We find that 



n 



(0) ^ ^(1) 







.(0) ^ ^(1) 



.(1) 



A'^-'^ ^ Jk) ^ (k+l) ^ ^ ^ ^ ^ (N) 



^0 ' 



(6.55) 



since (Tq is just the identity matrix. Note that if we had considered instead the product 

{af (6.56) 



then we would have obtained the same result. Therefore each term in P'^'^^ commutes with each of the 
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terms in P('). It follows that P^'^^ commutes with P('), and hence we can perform the decomposition 
of Eq. ( 6.50| ) in this case. 

we see 



Now consider the potential for the ss' — > pp' process when |m^ | = 1/2. Using Table 5.9 



that 



p(fc) 



pi'' 



P32 1^, 



(0) 

3 



(fc) 
'P33 



/ 



(0) 

3 



(6.57) 



and similarly for P^^\ Examining the product 



(6.58) 



we see that 



n 



JO 



(fc+1) 



(6.59) 



while if we had multiplied the factors in the opposite order we would have obtained 



(fe+i) 



J'-i) 



JO 



J'+i) 



JJV) 



(6.60) 



Since a^^a^"^ — *'''3°\ while ctj^VJ*''' = —iaf \ we see that the commutator of ^a^i^ ® '^3^'^ ^'^'^ 
|cr2°'' ® cj''! is not a c-number. Because the set of Ipq'^^I is not the same as the set of it 
follows that the commutator of P'^^^ and P'^'-' cannot be a c-number. Hence we cannot perform the 
separation of Eq. ( |6.5C| ) in this case. 



Considering the potential for the sp ps process and consulting Table 3.£, we see that the 
coefRcients of cri ® CT3 and 172 ® (Ti are both nonzero. Therefore the argument of the previous 



paragraph will hold here too, and again we cannot perform the separation of Eq. (6.50) 
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6.3.5 Averaging the ss' ^ pp' process when \mj\ = 3/2 for the p' atom 

In this case we can perform the separation of Eq. ( 6.5(]| ), and so we see from Eq. ( 6.44 ) that the 
averaging problem reduces the computation of the average of 



N 



2N- 



k=l 



-Tr n exp -/3^f pC^) 



Using Table 6.S, we can write 



(6.61) 



exp [-/JVfe^pC^-)] = exp 








= exp 




exp 


3 
b=l 



b=l 



(6.62) 



Ictq"^ ® '^o'^'' I i^^^^ 2^+-'^ X 2^+^ identity matrix, which certainly commutes with ^a^^^ ® cr^'^'' | . 



Here v^^p depends on rfc, and specifically 



Jk) _ A^A^ Jk) 

' k 



where p^^^ is the same quantity as in Eq. ( 6.53 ) . 
If we define 



fc=i 



and 



b=l 



then we can rewrite Eq. ( 6.62| ) as 



(6.63) 



(6.64) 



(6.65) 



exp 



exp 



exp I 



({af ® [-/JvC^) • aC^)] }) . (6.66) 



Now we note that 



oo ^ 

exp ((73 ® ^) = E ~j ('^a ® 

OO ^ 



3=0 
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°° ^ / /1^2fc , /i\2! + l 

gp (E> (A) ^ erg (g) (A) 

^ (2A:)l ^ 



fe=0 ' ' /=0 

(jQ (Xi cosh A + (Ts <E) sinh A, 



(2Z + 1) 



so that 



exp ({4°) ® [-/JvC^) . a^'^)] }) 



= ai?^ (g) cosh 



(73°'' (g) sinh 



where to save space we define 



A derivation very similar to that of Eq. (6.67) shows that 



and 



cosh (ctq ^ A) = (Tq ^ cosh A, 



sinh ((To ® A) = do (8) sinh A, 



and by repeated appHcation of these last two relations we finally arrive at 



{{4 



exp I <; ai°^ 



}) - {'S 



(0) 



) cosh 



1 173°'' (Xi sinh 



From page 166 of Ref. [B3| we have the familiar result that 



1 if n is even, 



From this it is easy to show that 



cosh -/^v^*^) • (T^'^) = cosh /SuC') , 
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and 



Thus we see that 



sinh 



^ -v(fe) . CT^*^) smh 



(k) 



(6.75) 



({' 



exp I <: 



}) 



cosh 
sinh 



(6.76) 



and hence 



exp 



-/3lf p(fe) 



exp 
exp 



cosh 
sinh 



(jW}. (6.77) 



We now note that because Tr ctq = if a is equal to 1, 2, or 3, it follows that when we expand 
the product 



and take the trace, only the term 



9A'+l„('=)„(0 



(6.79) 



survives. This together with Eqs. ( |6.61 ) and ( |6.77 ) implies that 



AT 



N 



2N+1 



Tr Yl exp -/JVfc^pC^) =11^^? 



cosh 



f3v 



k=l 



k=l 



We can now change variables so that 



N 



(k) 



cosh 



(3v 



(fe) 



\fc=i 



is equal to a product of N copies of 



(6.80) 



(6.81) 



drk rl I dOk sin Ok J cL0k exp 



cosh 



(k) 



(6.82) 



Chapter 6: The Effect of Spin 



115 



This integral, in turn, is equal to the sum of the integrals 



— J drk rl J dOk smOk J dc^k exp |-/3 \v'^^ - t^^'^'^J | , 



(6.83) 



and 



^" Jo ^0 



dru rl I dOi. 



k sinOk / dc^k exp {-/? + u^'^)] } 



Once again using the coefficients of Table 6.S, we find that 



(6.84) 



«6o' =^T^[7-3cos(20fc)], 



(6.85) 



and 



rl 3V3 



V5-3cos(26ifc). 



(6.86) 



The integrals of Eqs. (6.83) and (6.84) are now of a familiar type, and they can be done using 



the methods explained in Sections |2.3| and 2.5. Unfortunately the angular integrals have to be done 
numerically, but the result is that|^ 



^ / dru rl I dOk sin 9k 



27r 



bk exp 



cosh 



= 1 - 5.18061-^v/^. (6.87) 



It follows that 



where 



2JV+1 



N 

-Tr [| exp [-/3yfe2p(fc)j ^ (^^y^^,^^^^ 



k=X 



Vs-pm = 5.18061 /i/l' 



N 



(6.88) 



(6.89) 



Recalling Eq. ( 2.12 ), we see that when we used the simple fj.fi' /r^ interaction potential of Chap- 
ter |2[ we obtained 

4^3/2 



N 



-fifi' — = 7.42444 fifi'— 



(6.90) 



Hence we arrive at 



0.697778 u. 



(6.91) 



-"^The astute reader will be relieved to know that the quantities + f '''' and ^i^q ^ — ' '^''^ are both positive for all 
values of and 0^, and so there is no difficulty in applying the methods of Sections 2.3 and p.5| here. 
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Inserting our averaged result back into Eq. (6.44), we see that 



_Tr5(a))^ — - 



2N+ 



2a 2a 



d/3 exp \-(3 (a^ + A^) - 2vspin\/f3 



(6.92) 



Comparing with Eq. ( ^.27[ ), one immediately sees that Eq. (6.92) gives the same results as the signal 
computed in Chapter]^ but with the quantity Wspin appearing in place of v. 



6.4 Summary 

In this chapter we have considered the detailed effects of spin on the system we are studying. We 
were only able to fully average the ss' — > pp' process when \mj\ —3/2 for the p' atom, but we 



can see from Fig. 1.4 that the case where \mj\ — 1/2 for the p' atom is not very different. In the 
case that we could treat fully, we found that the effect of the spin was to change the value of the 
parameter w by a factor of 0.697778 from the simple case where we neglect the effects of spin and 
consider an angular independent 1/r'^ interaction potential. 

The consideration of the effects of spin becomes extremely important when a magnetic field is 
applied to the system. It was reported by Renn et al. that the introduction of a magnetic field 
greatly proliferates the number of nondegenerate states present in a room temperature Rydberg 
gas. One can imagine applying a magnetic field that is parallel to the detuning electric field and of 
sufficient strength to lift the rrij degeneracy. The form of the ss' — > pp' interaction would then be 



given by the appropriate matrix element of Table 6.2 or |6.3|, and the form of the sp ps interaction 



would similarly be found in Table 3.4, Thus we see that this is one possible scheme for selecting the 



angular form of the interaction potential. 



Chapter 7 



Numerical Simulations 



7.1 Introduction 

In this chapter we discuss numerical simulations of the system we are considering. We start by 
describing in detail how the simulations are carried out, then we present the results, including 
comparisons of the simulations with the experimental data. 

7.2 Review of theoretical model 

Before discussing exactly how the simulations are carried out, we review the basic elements of the 



theory presented at the beginning of Sections 2.2 and 4.2.1 of this thesis. We consider one atom at 



the origin, initially in the state s', in interaction with a gas of s atoms through a resonant ss' — > pp' 
dipole-dipole interaction V . We also allow for interaction among the unprimed atoms via an sp — > ps 



process mediated by another dipole-dipole potential U. This system is described by Eqs. (4.1a) and 



N 



(4.1b), which are 

i ciQ = Aao + ^VfeCfe, (7.1a) 
fe=i 

N 

i Ck = VkaQ + ^UkiCi. (7-lb) 

1=1 
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Here oq (A, t) is the amplitude of the state in which the atom at the origin is in state s' and aU 
other atoms are in state s, while (A, t) is the amplitude of the state in which the atom at the 
origin is in state p' and the atom at is in state p, while all the others remain in state s. The 
quantities Vk and Uki are the interaction potentials and A = ep> + ep — e^' — is the detuning from 
resonance. 

The simplest possibilities for the dipole-dipole interactions V and U are 

Vk = ~^, (7.2) 

and 

Uki = -^. (7.3) 

' kl 

Here fi is the dipole matrix element connecting the s andp states, and /i' is the dipole matrix element 
connecting the s' and p' states. We also define Uu = 0, so that the restriction I k is automatic in 



sums of the type appearing in Eq. (7. It). The atoms are assumed to be sufficiently cold that during 
the time scale of interest they move only a very small fraction of their separation, and therefore Vk 
and Uki can be taken to be independent of time. 



As was stated in Section 2.5.1, the induced dipoles in the magneto-optical trap become aligned 
if the applied electric field is sufficiently large that the Stark effect has passed from the quadratic to 
the linear regime. This is the case considered, for example, in Refs. and js^. For the numer- 
ical simulation results of this thesis we use the slightly more complicated dipole-dipole interaction 
potential corresponding to this case, which is 

Vk = -^P2{cosek), (7.4) 

where 9k is the angle between and the applied electric field. Similarly, Uki is of the form 

Uki = -^P2icoseki), (7.5) 

for I 7^ fc, and Ukk = 0. These formulas for V and U are also valid, in some cases, in the quadratic 
Stark regime if a magnetic field parallel to the electric field is used to align the total angular 
momentum of the atoms, as we now discuss more fully. 
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We will show below that, using interactions such as Eqs. (7.2) and (7.3), or Eqs. (7.4) and 
( [7.5[ ), we are able to reduce the computation of any quantity to the diagonalization of an {N + 1) x 
{N + 1) matrix. On the other hand, we saw in Chapter ^ that in the quadratic Stark regime 
the actual interaction potentials depend on the quantum numbers rrij of the atoms involved. A 



complete consideration of the spins would require the diagonalization of a 2 



N+l 



X 2 



N+l 



matrix in 



this case, which is computationally prohibitive to say the least. The problem can be reduced again to 
O {N X N) by applying a magnetic field that is parallel to the electric field and of sufficient strength 
to lift the mj degeneracy. The V interaction is then given by the appropriate matrix element of 



Table |6.2| or 6.2 , and the U interaction is similarly found in Table pA\ We see from these tables that 



for many matrix elements the angular dependence is the same as in Eqs. (7.4) and (7.5), which can 



therefore be taken as representative. We have also shown in Section 6.3.5 that, when a p' state with 
\mj\ 3/2 is involved, each s atom interacts separately with the s' atom and, as the system evolves, 
the average number of p' states is the same as for a simple spin-independent interaction, but with 



an effective Vk- The resulting effective interaction energy Uspin is given by Eq. (3.91) to be 0.698 u, 



where v is the effective interaction energy for the interaction potential of Eq. (7.2) and is defined 
expHcitly in Eq. ( |2.12| ). For comparison, the effective interaction energy for the interaction potential 
of Eq. (7.4) is given by Eq. ( 2.75| ) to be 0.770 Thus we see that the exact angular form of the 
interaction potential for V does not greatly affect the results. The numerical simulations reported 
in this chapter show that the same is true for the angular form of U. 



7.3 Description of method 

If we consider a column vector C (A, t) = (ao (A, t) , ci (A, t) , C2 (A, t) , . . . ,cn (A, t))'^ , then clearly 



we can construct a matrix H such that the set of iV + 1 coupled differential equations of Eqs. (7.1a) 



and (7.1b) can be written in the form i C ~ HC. If we let {A„j} and {9m} be the eigenvalues and 
eigenvectors (represented as column vectors) of the real Hermitian matrix H, respectively, and if we 
further define e„ to be the column vector of length + 1 with the nth entry set equal to one and 
all others set equal to zero, then we can write e„ = X]m=0'^™"^™ ^^"^ ^™ ~ 12n=o (^^"1-^^ ' where 
Pmn = ()fn^n- The systcm starts off in the state described by cq, so it follows that 



C(A,t) 



= exp (—iHt) Co 



Chapter 7: Numerical Simulations 



120 



N 



= ^ PmO exp {-iX,nt)0m 
m=0 
N N 

= X! X! l^"^oPnm exp {-iXmt) e„. (7.6) 

m— n— 

The quantity oq (A,<) is just the cq component of this expression, so 

N 

ao (A, i) = eJC (A, t) = ^1 /^"o/^o™ ^xp {^iX^t) . (7.7) 

m=0 

Thus we see that, in this form, solving for ao at a given time becomes just a matter of diagonahzing 
the matrix H. 

To carry out the numerical simulations, we first generate a set of random positions for the 



unprimed atoms. For the results presented in this Sections 7.5 and 7.6, these random positions are 
uniformly distributed so that the results can be more easily compared with the analytical results 
of previous chapters. However, when proceeding numerically it is just as easy to compute random 
positions for a gas blob with any given density profile. For the experiments we are modeling ^, 
the density profile is a Gaussian with an aspect ratio of approximately 1:1:5 |23|, and so this is 



what we use to generate the results of Section 7.7. Once the atomic positions are generated, they 



are used to compute the interaction potentials Vk and Uki, and these are in turn used to compute 
the matrix H. We then diagonalize H and construct the quantity ao at a set of values of t using 



Eq. (7.7). We do this for the same set of times for many random distributions of the imprimed atoms 



and thereby obtain an average of the quantities ao (A, t) and |ao (A, t)|^ at a given set of times over 
many configurations of the system. The signal measured by the experiments is then proportional to 
the average of 1- |ao {A,t)\. We can further obtain plots of the lincwidth as a function of time by 
generating signal versus time curves for many values of the detuning. A, and then determining for 
each value of time which of these signals is equal to one half the value of the resonance signal. One 



can also directly generate Fourier transforms from Eq. (7.7), as we did in Section 4.6.2. 

By comparing the results of the simulations to the analytical results of Chapter ^ we find that 
an average over 10^ realizations, each consisting of a single s' atom and 100 s atoms, is sufficient, as 



can be seen in Fig. 7.1. These are the parameters we have used for all the simulation plots presented 
in this thesis, although wc have run the simulations for as many as 800 s atoms. 

The simulations are written in Fortran 90, and the computation process is greatly accelerated 
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by the use of a Beowulf cluster of thirteen Linux machines. Using this chister, the computation can 
be carried out in parallel, as we now describe. A master program is started on one of the machines. 
The master then starts up a slave program on each machine in the cluster. The master initializes 
the slaves by sending them seeds for their random number generators.]^ The slaves then generate 
random distributions for the atoms, compute whatever quantity is required, and send the result 
back to the master. The master collects and averages the results returned by the slaves and saves 
the result to a data file. The master also takes care of any interaction with the user. The passing of 
messages and data between the master and slave processes is handled through the software program 
Parallel Virtual Machine (PVM) 

We use the (parallel) random number generation routines gasdev and rani described in Refs. ||4^ 
and to generate the random distributions of atoms. The routine gasdev generates Gaussian 
deviates with unit standard deviation, while rani generates uniform deviates between zero and one. 
The computation of eigenvalues and eigenvectors for the matrix H corresponding to each realization 
is performed using the ssyevx routine of the LAPACK package, which is described in Ref. |^. 

7.4 The effect of many s' atoms 
7.4.1 The chemical argument 

Considering the ss' pp' process as a chemical reaction, we have the rate equation 

^ = k{NsN,,-NpNp,), (7.8) 

where /c is a constant with the units of inverse time, which is equivalent to energy since we take h 
equal to unity. There are similar rate equations for Np, Ng, and N^', but it is sufficient to consider 
only one of these because of the following constraints. 

We know that Np — Npi , since if there is a p or p' atom present it can only have been produced 
as a result of a pair of atoms in the s and s' states making the ss' —* pp' transition. We also know 
that N = Ns + Np ^ Ns + Npi and N' — Ng' + Npi , since primed and unprimed atoms can change 

^It is vital that the master control the random number generator seeds, since each slave program must be seeded 
differently so that each generates a different set of realizations for the atomic positions. 
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from s to p states and vice versa, but can never be created or destroyed. Therefore we have 



dNp, 
dt 



= k[iN~Np){N' -Np,)-NpNp,] 
= k[iN~Np,)iN' -Np,)~Np,Np- 
^ k[NN' -{N + N')Np>]. 



(7.9) 



With the initial condition Npi = at i = 0, this gives 



Np' = 



NN' 
N + N' 



1 - e 



-[N+N')kt 



and at small times 



7V„, = NN'kt. 



(7.10) 



(7.11) 



7.4.2 Extrapolation to many s' atoms 

In our case, the reaction rate is not constant and the approach to equilibrium is not accurately 



described by Eq. ( 7.1C ). However, the result we have from Chapter g in the limit N' <^ N starts off 
at small times as 

(7.12) 



Np, = c^i^i'—N't, 



where C is a numerical constant that does not concern us here. This expression is fine for all values 
of N and N\ and in particular satisfies the requirement that Np — Np'. It is of the same form as 



Eq. (7.11) and shows that at small times we have 



k = C 



(7.13) 



Furthermore, we can write 



NN' / ,N + N' 

Np, = F ijLii' 1, . 

P N + N' 



n 



(7.14) 



where the dots indicate that the function F depends on other dimensionless ratios that character- 
ize the system but do not involve t. The prefactor NN' / {N + N') and the characteristic energy 



fifi' {N + N') /Q. are suggested by the chemical argument of Section 7.4.1, since these factors are 
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explicitly seen in Eq. (7.10) if we keep in mind Eq. ( 7.13| ). The full function F {x, . . .) is computable 



only in the dilute limits N' ^ N oi N ^ N' , but we know in general that it starts off linearly in x. 
For large values of t we have, again for N' ^ N, 

where / (y, z) is some function that can be computed numerically using the methods of Section [7.3| . 
This result cannot be correct for all values of N/N' . In fact, for N <^ N' we can compute Np and 
we find 

Np = Nf(^,^]. (7.16) 
Since Npr is equal to Np, we see that this (correct) result is quite different from what one would 



obtain using Eq. ( 7.15 ), which is incorrect in this limit. Interpolating between the limits given by 



Eqs. ( 7.15 ) and ( 7.16 ), we postulate that for large t 



We can join this result smoothly to the small t result of Eq. ( 7.12 ) by postulating a general result 
of the form 

NN' ( ,N + N' A N ^Ji iV' 

where the fimction F{x,y,z) is proportional to x for small x and approaches f{y,z) for large x. 
We have also used the same scaling for the quantities t and A, so that xy — At. 

With this approach we can use the function F{x,y,z), which can be computed numerically in 
the dilute limit, as a reasonable guess for the functional dependence on the suitable dimensionless 
variables for any value of N' /N. 

The combination 

N II N' a' , , 

^ • ' (7.19) 



N + N' N + N' fi 

is just an interpolation between the correct results in the limits of small N and small N' . It is 
reasonable, however, since the quantities 

iV/x^f, (7.20) 
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and 



N'fifi'- (7.22) 



give the effective U for the sp — > ps interaction for a fraction N' / {N + N') of the atoms, the effective 
U' for the s'p' p's' interaction for a fraction N/ {N + N') of the atoms, and the V interaction for 



the ss' — > pp' process for ah the atoms. The quantity of Eq. (7.19) is obtained by taking the ratio 

\, " ■ 7.23 

7.5 Results of simulations for the signal 



For the results of this section, we well as those of Section 7.6, we use the interaction potentials of 



Eqs. (7A) and (7^). We also use a uniform distribution to generate the random atomic positions, 
so that the results can be more easily compared with the analytical results of previous chapters. 

Figs. [7.l| - |7.5| display the numerical simulations of 1 — ^|ao|^^, corresponding to the experimental 
signal Np' /N' , as a function of vt for u/v ~ 0, 1/4, 1, 2, and 4, and enough values of the detuning A 
to give an idea of the resonance width in each case. The FWHM at large t can be directly deduced 
from the intercepts of the graphs with the right edge of the figure, as is discussed more fully below. 

In addition. Fig. shows the exact results for u — that were discussed in Chapter Com- 
parison with the numerical simulations, which were carried out for N = 100, shows that the infinite 
system is already well simulated, but that noticeable wiggles remain after averaging over 10^ realiza- 
tions. The reason for this persistence of fluctuations is that most realizations do not look at all like 
the average. In fact, the variance with realization is comparable to the average. In particular, for 
small t the average is given by (y^/2) vt while each realization is quadratic in t, as already discussed 
in Ref. It can be seen from the graphs that the linear rise law, {y/TT/2) vt, is independent both 
of A and of u and holds for a substantial range of vt. 
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Figure 7.1; Plots of 1 — (^\ao (A,t)|^^ for u = and several values of the detuning A. The exact results of Chapter 
are shown for comparison. 
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Figure 7.2: Plots of 1 — (\ao (A,t)|^\ for u = 0.25v and several values of the detuning A. 
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Figure 7.3: Plots of 1 — (\ao (A,t)|^\ for u = v and several values of the detuning A. 
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Figure 7.4: Plots of 1 - 



(|ao (A,t)|2^ foru 



= 2v and several values of the detuning A. 




Figure 7.5: Plots of 1 — (\ao (A,t)|^\ for u = 4v and several values of the detuning A. 
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7.6 Results of simulations for the width 



The time-dependent resonance lineshape is obtained from the ordinates of the curves in Figs. 7.1-7.5 
at a given vt. In general the Hneshape is not Lorentzian, but we never find a spHt hue with a minimum 
at resonance as reported in Rcf. |4^. Thus it is fair to characterize the hne by its FWHM. This 
cannot be directly compared with an experimental width, which is usually the FWHM of a Lorentzian 
fitted to a noisy experimental line, but it is indicative of trends. For a more accurate comparison 
with experiment, one should of course convolute the calculated lineshape with the broadening due 
to other effects. 



Fig. 7.6 shows one half of the FWHM, w, as a function of vt for u — iv. One sees the typical 
decrease of w from a universal t~^ behavior at small vt towards an asymptotic value, which is almost 
reached at vt = 15. The t^^ law, which results simply from the uncertainty principle H, has 
already been discussed in Ref. where on the basis of a small A expansion it was estimated that 
w = \/l2/< for small t. This t^^ behavior was observed experimentally by Renn et al. [ |60[ in a 
different system, and was later observed by Anderson et al. |j in the Rb system. 



The near-asymptotic value of the FWHW at vt — 15 is shown in Fig. 7.7 as a function of u/v. 
Overall, the graph is roughly linear except at small u/v, although a careful examination shows that 
there are really two linear regimes with somewhat different slopes that meet at u/v w 3.5. There is 
an unexpected minimum around u = v /A, which is not a numerical artifact, but can be seen directly 



from Figs. 7.1-7.5 by comparing the plot for u — w/4 with those for m = and u = v. What happens 
is that the resonance line always grows taller and wider as u increases, but for small u its height 
(the signal at A = 0) grows faster than its width. 

7.7 Comparison with experimental data 



For the results in this section, we use the interaction potentials of Eqs. ( |7.4| ) and (7.5). Because 
it corresponds more closely to the experimental situation, we also use a Gaussian distribution with 
aspect ratio 1 : 1 : 5 to generate the random atomic positions. 



Figs. 7.S and 7.9 show the experimental signal as a function of time, together with a signal 
computed using the numerical simulations and fit to the data. The experimental data for both 
graphs are the data of Fig. 3.8 on page 71 of Ref. [Q, but we show here all the data points. 

The fits are made by first determining the vertical coordinates where the experimental signal 
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Figure 7.6: A plot of the width as a function of time for u = 4v. 
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Figure 7.7: A plot of the width at lit = 15 as a function of the ratio u/v. 
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starts and stops. The difference between the two gives the vertical scale of the experimental signal, 
which is necessary because as one can see the experimental signal does not simply start at zero 



and rise to some value less than one. For example, in Fig. 7.i we assumed that the experimental 



signal started at 900 and ended at 1850 so that the vertical scale was 950. The initial slope is then 



determined by fitting the first several points to a straight line. For Fig. 7.8 we obtained an initial 
slope of 4107 fj,s~^. We then have 

= 3.83115 /is-i, (7.24) 

which gives a value for the parameter v and hence the density N/il. This fixes the horizontal scale, 
and so all that remains is to adjust the parameter u/v until we find a good fit. 



Fig. 7.1C shows the experimental width as a function of time. The data of Anderson et al. are 
taken from Fig. 4-18 on page 149 of Ref. while the data of Lowell et al. are taken from Fig. 3.10 
on page 75 of Ref |4^ . The discrepancy between the two sets of data has not yet been resolved by 
the experimentalists. Although the result of the simulations fit Anderson's data quite well for the 
parameters u — and N/Vl ~ 2 y. 10^ cm^"^, it is important to note that these data are for rather 
small times, and are within or very close to the region where the universal t^^ behavior discussed in 



Section 7.6 is present. Thus for a wide range of values of the parameter u there will exist densities 



N/n that produce a good fit to the Anderson data. 



7.8 Summary 



We have seen in this chapter how to develop numerical simulations to solve Eqs. ( 7.1a ) and ( 7.1b ) 
for any quantity we want and average the result over many realizations of the system. Although 
Eqs. (7.1a) and ( 7.1b| ) apply to the sparse case where we begin initially with a single s' atom in a 



sea of s atoms, we have also discussed in this chapter how to extend these results to the nonsparse 
case. The simulation results agree very well with the exact analytical results of Chapter |^, as well 
as with the approximate analytical results of Chapter ^ The results of the simulations can also be 
compared with the experimental data of Anderson et al. ^ and Lowell et al. pO[ , and we find 
that the agreement is again quite good. 
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Figure 7.9: The experimental signal as a function of time (triangles), together with a signal curve computed using 
the numerical simulations and fit to the data. The fit parameters are u = 12v and N/Q = 6.70 x 10* cm~*. 
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Figure 7.10: The experimental width as a function of time. The circles are the data of Anderson et al. g|, while 
the squares are the data of Lowell po| ]. The triangles are a width curve computed using the numerical simulations 
and fit to the data. The fit parameters are u = and N/Q, = 2 X 10® cm~^. 
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Final Thoughts 

Although we have focused chiefly on resonant dipole-dipole interactions among frozen Rydberg atoms 
in this thesis, there are many connections of this work with related experiments, as well as more 
general applications of the techniques we have developed. 

In this thesis we have always considered the atoms to be frozen in place. It is true that the atoms 
are sufficiently cold that they move only a small fraction of their average interatomic spacing during 
the course of the experiment, but close pairs of atoms can be separated by distances considerably 
smaller than this, and hence the effects of interatomic motion can be considered. Certainly this would 
be a worthwhile endeavor, since one would ideally like to eventually have a theory that reduces to 
the frozen results of this thesis at low temperatures and, as the temperature increases, eventually 
makes the transition to the high-temperature regime where the results are dominated by binary 
atomic collisions. Some ideas for the inclusion of finite temperature effects can be gleaned from the 
extensive literature on dipolar liquids jl^, |2^, |3^, 

Dipole-dipole interactions among frozen Rydberg atoms may one day be useful for quantum 
computing, as methods for producing entanglement in such systems have already been studied 
theoretically |3^, ^ as well as experimentally This method for quantum computing has 

advantages over others that have been proposed. Entanglement of photon pairs, for example, has 
the disadvantage that the production of such pairs is not very efficient. 

Experiments have been done recently that illustrate the creation of ultracold plasmas using 
Rydberg atoms jsj, |3^, and even the spontaneous evolution of Rydberg atoms into an ultracold 
plasma |52| . Unlike in a high-temperature plasma where the kinetic energy of the particles always 
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dominates, an ultracold plasma can be strongly coupled because the Coulomb interaction between 
neighboring particles can exceed their kinetic energy. Therefore such plasmas are of theoretical in- 
terest because they exhibit qualitatively different features than their high-temperature counterparts. 
Although we have averaged over random distributions of atoms in this work, it is also possible 



to construct periodic optical potentials 29, ^ |61|| and even quasiperiodic optical lattices [ p7[ . 
In the periodic case it is necessary to, instead of averaging, consider instead the appropriate dipolar 
sums for regular lattices, which are discussed in Refs. jl^ and ||20|]. 
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